[Pomp-commits] r490 - in pkg: . R data inst inst/doc man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri May 20 14:19:27 CEST 2011
Author: kingaa
Date: 2011-05-20 14:19:27 +0200 (Fri, 20 May 2011)
New Revision: 490
Modified:
pkg/DESCRIPTION
pkg/R/mif.R
pkg/R/pfilter.R
pkg/R/pmcmc.R
pkg/data/blowflies.rda
pkg/data/dacca.rda
pkg/data/euler.sir.rda
pkg/data/gillespie.sir.rda
pkg/data/gompertz.rda
pkg/data/ricker.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/ChangeLog
pkg/inst/NEWS
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/gompertz-multi-mif.rda
pkg/inst/doc/gompertz-trajmatch.rda
pkg/inst/doc/intro_to_pomp.pdf
pkg/inst/doc/ricker-mif.rda
pkg/inst/doc/ricker-probe-match.rda
pkg/man/pfilter.Rd
pkg/src/pfilter.c
pkg/tests/ou2-forecast.R
pkg/tests/ou2-forecast.Rout.save
pkg/tests/ou2-mif.R
pkg/tests/ou2-mif.Rout.save
pkg/tests/ou2-procmeas.R
pkg/tests/ou2-procmeas.Rout.save
pkg/tests/ou2-trajmatch.Rout.save
pkg/tests/pfilter.R
pkg/tests/pfilter.Rout.save
pkg/tests/sir-icfit.R
pkg/tests/sir-icfit.Rout.save
pkg/tests/synlik.R
pkg/tests/synlik.Rout.save
Log:
- pfilter has been redesigned to allow variable numbers of particles
- pmcmc and mif have been modified to accommodate this
- the 'saved.states' and 'saved.params' slots in 'pfilterd.pomp' objects are now lists of rank-2 arrays (they were each a rank-3 array)
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/DESCRIPTION 2011-05-20 12:19:27 UTC (rev 490)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.37-2
-Date: 2011-05-19
+Version: 0.38-1
+Date: 2011-05-22
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing, Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
URL: http://pomp.r-forge.r-project.org
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/R/mif.R 2011-05-20 12:19:27 UTC (rev 490)
@@ -107,12 +107,27 @@
if (missing(particles))
stop("mif error: ",sQuote("particles")," must be specified",call.=FALSE)
-
+
+ ntimes <- length(time(object))
if (missing(Np))
stop("mif error: ",sQuote("Np")," must be specified",call.=FALSE)
+ if (is.function(Np)) {
+ Np <- try(
+ vapply(seq.int(from=0,to=ntimes,by=1),Np,numeric(1)),
+ silent=FALSE
+ )
+ if (inherits(Np,"try-error"))
+ stop("if ",sQuote("Np")," is a function, it must return a single positive integer")
+ }
+ if (length(Np)==1)
+ Np <- rep(Np,times=ntimes+1)
+ else if (length(Np)!=(ntimes+1))
+ stop(sQuote("Np")," must have length 1 or length ",ntimes+1)
+ if (any(Np<=0))
+ stop("number of particles, ",sQuote("Np"),", must always be positive")
+ if (!is.numeric(Np))
+ stop(sQuote("Np")," must be a number, a vector of numbers, or a function")
Np <- as.integer(Np)
- if ((length(Np)!=1)||(Np < 1))
- stop("mif error: ",sQuote("Np")," must be a positive integer",call.=FALSE)
if (missing(ic.lag))
stop("mif error: ",sQuote("ic.lag")," must be specified",call.=FALSE)
@@ -136,20 +151,6 @@
if (Nmif<0)
stop("mif error: ",sQuote("Nmif")," must be a positive integer",call.=FALSE)
- if (verbose) {
- cat("performing",Nmif,"MIF iteration(s) to estimate parameter(s)",
- paste(pars,collapse=", "))
- if (length(ivps)>0)
- cat(" and IVP(s)",paste(ivps,collapse=", "))
- cat(" using random-walk with SD\n")
- print(rw.sd)
- cat(
- "using",Np,"particles, variance factor",var.factor,
- "\ninitial condition smoothing lag",ic.lag,
- "and cooling factor",cooling.factor,"\n"
- )
- }
-
theta <- start
sigma <- rep(0,length(start))
@@ -203,7 +204,7 @@
## initialize the particles' parameter portion...
P <- try(
- particles(tmp.mif,Np=Np,center=theta,sd=sigma.n*var.factor),
+ particles(tmp.mif,Np=Np[1],center=theta,sd=sigma.n*var.factor),
silent=FALSE
)
if (inherits(P,'try-error'))
Modified: pkg/R/pfilter.R
===================================================================
--- pkg/R/pfilter.R 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/R/pfilter.R 2011-05-20 12:19:27 UTC (rev 490)
@@ -9,8 +9,8 @@
filter.mean="array",
eff.sample.size="numeric",
cond.loglik="numeric",
- saved.states="array",
- saved.params="array",
+ saved.states="list",
+ saved.params="list",
seed="integer",
Np="integer",
tol="numeric",
@@ -39,22 +39,40 @@
if (length(params)==0)
stop(sQuote("pfilter")," error: ",sQuote("params")," must be specified",call.=FALSE)
- if (missing(Np))
- Np <- NCOL(params)
-
if (missing(tol))
stop(sQuote("pfilter")," error: ",sQuote("tol")," must be specified",call.=FALSE)
one.par <- FALSE
times <- time(object,t0=TRUE)
ntimes <- length(times)-1
+
+ if (missing(Np))
+ Np <- NCOL(params)
+ if (is.function(Np)) {
+ Np <- try(
+ vapply(seq.int(from=0,to=ntimes,by=1),Np,numeric(1)),
+ silent=FALSE
+ )
+ if (inherits(Np,"try-error"))
+ stop("if ",sQuote("Np")," is a function, it must return a single positive integer")
+ }
+ if (length(Np)==1)
+ Np <- rep(Np,times=ntimes+1)
+ else if (length(Np)!=(ntimes+1))
+ stop(sQuote("Np")," must have length 1 or length ",ntimes+1)
+ if (any(Np<=0))
+ stop("number of particles, ",sQuote("Np"),", must always be positive")
+ if (!is.numeric(Np))
+ stop(sQuote("Np")," must be a number, a vector of numbers, or a function")
+ Np <- as.integer(Np)
+
if (is.null(dim(params))) {
one.par <- TRUE # there is only one parameter vector
coef(object,names(params)) <- unname(params) # set params slot to the parameters
params <- matrix(
params,
nrow=length(params),
- ncol=Np,
+ ncol=Np[1],
dimnames=list(
names(params),
NULL
@@ -71,21 +89,13 @@
## set up storage for saving samples from filtering distributions
if (save.states)
- xparticles <- array(
- data=NA,
- dim=c(nvars,Np,ntimes),
- dimnames=list(statenames,NULL,NULL)
- )
+ xparticles <- vector(mode="list",length=ntimes)
else
- xparticles <- array(dim=c(0,0,0))
+ xparticles <- list()
if (save.params)
- pparticles <- array(
- data=NA,
- dim=c(length(paramnames),Np,ntimes),
- dimnames=list(paramnames,NULL,NULL)
- )
+ pparticles <- vector(mode="list",length=ntimes)
else
- pparticles <- array(dim=c(0,0,0))
+ pparticles <- list()
random.walk <- !missing(.rw.sd)
if (random.walk) {
@@ -149,12 +159,12 @@
filt.m <- array(dim=c(0,0))
for (nt in seq_len(ntimes)) {
-
+
## advance the state variables according to the process model
X <- try(
rprocess(
object,
- x=x,
+ xstart=x,
times=times[c(nt,nt+1)],
params=params,
offset=1
@@ -208,7 +218,7 @@
xx <- try(
.Call(
pfilter_computations,
- X,params,
+ X,params,Np[nt+1],
random.walk,sigma,
pred.mean,pred.var,
filter.mean,one.par,
@@ -242,11 +252,11 @@
}
if (save.states) {
- xparticles[,,nt] <- x
+ xparticles[[nt]] <- x
}
if (save.params) {
- pparticles[,,nt] <- params
+ pparticles[[nt]] <- params
}
if (verbose && ((ntimes-nt)%%5==0))
Modified: pkg/R/pmcmc.R
===================================================================
--- pkg/R/pmcmc.R 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/R/pmcmc.R 2011-05-20 12:19:27 UTC (rev 490)
@@ -93,11 +93,26 @@
if (missing(dprior.fun))
stop("pmcmc error: ",sQuote("dprior")," must be specified",call.=FALSE)
+ ntimes <- length(time(object))
if (missing(Np))
stop("pmcmc error: ",sQuote("Np")," must be specified",call.=FALSE)
+ if (is.function(Np)) {
+ Np <- try(
+ vapply(seq.int(from=0,to=ntimes,by=1),Np,numeric(1)),
+ silent=FALSE
+ )
+ if (inherits(Np,"try-error"))
+ stop("if ",sQuote("Np")," is a function, it must return a single positive integer")
+ }
+ if (length(Np)==1)
+ Np <- rep(Np,times=ntimes+1)
+ else if (length(Np)!=(ntimes+1))
+ stop(sQuote("Np")," must have length 1 or length ",ntimes+1)
+ if (any(Np<=0))
+ stop("number of particles, ",sQuote("Np"),", must always be positive")
+ if (!is.numeric(Np))
+ stop(sQuote("Np")," must be a number, a vector of numbers, or a function")
Np <- as.integer(Np)
- if ((length(Np)!=1)||(Np < 1))
- stop("pmcmc error: ",sQuote("Np")," must be a positive integer",call.=FALSE)
if (missing(hyperparams))
stop("pmcmc error: ",sQuote("hyperparams")," must be specified",call.=FALSE)
Modified: pkg/data/blowflies.rda
===================================================================
(Binary files differ)
Modified: pkg/data/dacca.rda
===================================================================
(Binary files differ)
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gillespie.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/gompertz.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ricker.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 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/inst/ChangeLog 2011-05-20 12:19:27 UTC (rev 490)
@@ -1,5 +1,47 @@
+2011-05-19 kingaa
+
+ * [r489] tests/ou2-procmeas.R, tests/ou2-procmeas.Rout.save: - new
+ test for dmeasure and dprocess in ou2
+ * [r488] data/ou2.rda, inst/data-R/ou2.R, man/ou2.Rd, src/ou2.c,
+ tests/ou2-probe.Rout.save: - rework the OU2 model slightly
+ * [r487] tests/sir.Rout.save: - update sir test
+ * [r486] DESCRIPTION, data/blowflies.rda, data/dacca.rda,
+ data/gompertz.rda, data/ou2.rda, data/ricker.rda, data/rw2.rda,
+ data/verhulst.rda, inst/NEWS, inst/doc/gompertz-multi-mif.rda,
+ inst/doc/gompertz-trajmatch.rda, inst/doc/intro_to_pomp.Rnw,
+ inst/doc/intro_to_pomp.pdf, inst/doc/ricker-mif.rda,
+ inst/doc/ricker-probe-match.rda, tests/sir-icfit.Rout.save,
+ tests/sir.R: - minor improvements to 'intro to pomp' vignette
+ - update data objects
+ - update tests
+ * [r485] data/euler.sir.rda, data/gillespie.sir.rda,
+ inst/examples/sir.c, src/sir.c: - change the measurement model
+ for the SIR example
+ * [r484] inst/data-R/ou2.R, inst/doc/advanced_topics_in_pomp.Rnw,
+ inst/doc/advanced_topics_in_pomp.pdf: - improve the 'advanced
+ topics' vignette
+ - update the data()-loadable ou2 object
+ * [r483] inst/examples/ou2.c, src/ou2.c: - modifications to the C
+ codes for the OU2 example
+ * [r482] R/plugins.R: - change 'discrete.time.sim' so that
+ 'delta.t' is not a required argument of 'step.fun'
+ * [r481] man/pomp-package.Rd, man/pomp.Rd: - improve the
+ documentation a bit
+
+2011-05-17 kingaa
+
+ * [r480] man/pomp.Rd: - put in links to pomp methods
+
+2011-05-16 kingaa
+
+ * [r479] DESCRIPTION, R/sannbox.R, R/sobol.R, R/spect-match.R,
+ R/spect.R: - replace some instances of 'sapply' with 'vapply'
+
2011-05-13 kingaa
+ * [r478] DESCRIPTION, inst/ChangeLog,
+ inst/doc/advanced_topics_in_pomp.pdf, inst/doc/intro_to_pomp.pdf:
+ - update DESCRIPTION and ChangeLog
* [r477] data/blowflies.rda, data/dacca.rda, data/euler.sir.rda,
data/gillespie.sir.rda, data/ou2.rda, data/rw2.rda,
data/verhulst.rda, inst/NEWS,
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/inst/NEWS 2011-05-20 12:19:27 UTC (rev 490)
@@ -1,6 +1,11 @@
NEWS
0.37-2
+ o In 'pfilterd.pomp' objects, 'saved.states' and 'saved.params' slots are now length-ntimes lists of arrays.
+
+ o It is now possible to use variable numbers of particles in 'pfilter'.
+ The 'Np' argument may be specified as a single number as before, or now as either a vector of numbers or a function.
+
o The "advanced topics" vignette has been much augmented and improved.
o The 'euler.sir' example (load it with 'data(euler.sir)') has a different measurement model: a discretized normal distribution with the same mean and variance as before.
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/gompertz-multi-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/gompertz-trajmatch.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ricker-mif.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/ricker-probe-match.rda
===================================================================
(Binary files differ)
Modified: pkg/man/pfilter.Rd
===================================================================
--- pkg/man/pfilter.Rd 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/man/pfilter.Rd 2011-05-20 12:19:27 UTC (rev 490)
@@ -86,13 +86,14 @@
numeric vector containing the conditional log likelihoods at each time point.
}
\item{saved.states}{
- If \code{pfilter} was called with \code{save.states=TRUE}, this is the array of state-vectors at each time point, for each particle.
- An array with dimensions \code{nvars}-by-\code{Np}-by-\code{ntimes}.
- In particular, \code{saved.states[,i,t]} can be considered a sample from \eqn{f[X|y_{1:t}]}.
+ If \code{pfilter} was called with \code{save.states=TRUE}, this is the list of state-vectors at each time point, for each particle.
+ It is a length-\code{ntimes} list of \code{nvars}-by-\code{Np} arrays.
+ In particular, \code{saved.states[[t]][,i]} can be considered a sample from \eqn{f[X_t|y_{1:t}]}.
}
\item{saved.params}{
- If \code{pfilter} was called with \code{save.params=TRUE}, this is the array of parameter-vectors at each time point, for each particle.
- It is an array with dimensions \code{npars}-by-\code{Np}-by-\code{ntimes}.
+ If \code{pfilter} was called with \code{save.params=TRUE}, this is the list of parameter-vectors at each time point, for each particle.
+ It is a length-\code{ntimes} list of \code{npars}-by-\code{Np} arrays.
+ In particular, \code{saved.params[[t]][,i]} is the parameter portion of the i-th particle at time \eqn{t}.
}
\item{seed}{
the state of the random number generator at the time \code{pfilter} was called.
Modified: pkg/src/pfilter.c
===================================================================
--- pkg/src/pfilter.c 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/src/pfilter.c 2011-05-20 12:19:27 UTC (rev 490)
@@ -3,26 +3,26 @@
#include "pomp_internal.h"
#include <Rdefines.h>
-static void nosort_resamp (int n, double *w, int *p, int offset)
+static void nosort_resamp (int nw, double *w, int np, int *p, int offset)
{
int i, j;
double du, u;
- for (j = 1; j < n; j++) w[j] += w[j-1];
+ for (j = 1; j < nw; j++) w[j] += w[j-1];
- if (w[n-1] <= 0.0)
+ if (w[nw-1] <= 0.0)
error("non-positive sum of weights");
- du = w[n-1] / ((double) n);
+ du = w[nw-1] / ((double) np);
u = runif(-du,0);
- for (i = 0, j = 0; j < n; j++) {
+ for (i = 0, j = 0; j < np; j++) {
u += du;
while (u > w[i]) i++;
p[j] = i;
}
if (offset) // add offset if needed
- for (j = 0; j < n; j++) p[j] += offset;
+ for (j = 0; j < np; j++) p[j] += offset;
}
@@ -36,7 +36,7 @@
GetRNGstate();
n = LENGTH(weights);
PROTECT(perm = NEW_INTEGER(n)); nprotect++;
- nosort_resamp(n,REAL(weights),INTEGER(perm),1);
+ nosort_resamp(n,REAL(weights),n,INTEGER(perm),1);
UNPROTECT(nprotect);
PutRNGstate();
return(perm);
@@ -49,7 +49,7 @@
// weights are used in filtering mean computation.
// if length(weights) == 1, an unweighted average is computed.
// returns all of the above in a named list
-SEXP pfilter_computations (SEXP x, SEXP params,
+SEXP pfilter_computations (SEXP x, SEXP params, SEXP Np,
SEXP rw, SEXP rw_sd,
SEXP predmean, SEXP predvar,
SEXP filtmean, SEXP onepar,
@@ -62,9 +62,9 @@
SEXP retval, retvalnames;
double *xpm = 0, *xpv = 0, *xfm = 0, *xw = 0, *xx = 0, *xp = 0;
SEXP dimX, dimP, newdim, Xnames, Pnames, pindex;
- int *dim, *pidx, lv;
+ int *dim, *pidx, lv, np;
int nvars, npars = 0, nrw = 0, nreps, offset, nlost;
- int do_rw, do_pm, do_pv, do_fm, is_op, all_fail = 0;
+ int do_rw, do_pm, do_pv, do_fm, do_par_resamp, all_fail = 0;
double sum, sumsq, vsq, ws, w, toler;
int j, k;
@@ -81,14 +81,16 @@
error("'states' and 'params' do not agree in second dimension");
PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
+ np = INTEGER(AS_INTEGER(Np))[0]; // number of particles to resample
+
PROTECT(rw_names = GET_NAMES(rw_sd)); nprotect++; // names of parameters undergoing random walk
do_rw = *(LOGICAL(AS_LOGICAL(rw))); // do random walk in parameters?
do_pm = *(LOGICAL(AS_LOGICAL(predmean))); // calculate prediction means?
do_pv = *(LOGICAL(AS_LOGICAL(predvar))); // calculate prediction variances?
do_fm = *(LOGICAL(AS_LOGICAL(filtmean))); // calculate filtering means?
- is_op = *(LOGICAL(AS_LOGICAL(onepar))); // are all cols of 'params' the same?
- is_op = is_op && !do_rw;
+ do_par_resamp = *(LOGICAL(AS_LOGICAL(onepar))); // are all cols of 'params' the same?
+ do_par_resamp = !do_par_resamp || do_rw || (np != nreps); // should we do parameter resampling?
PROTECT(ess = NEW_NUMERIC(1)); nprotect++; // effective sample size
PROTECT(loglik = NEW_NUMERIC(1)); nprotect++; // log likelihood
@@ -217,19 +219,19 @@
if (!all_fail) { // resample the particles unless we have filtering failure
int xdim[2];
- int sample[nreps];
+ int sample[np];
double *ss, *st, *ps, *pt;
// create storage for new states
- xdim[0] = nvars; xdim[1] = nreps;
+ xdim[0] = nvars; xdim[1] = np;
PROTECT(newstates = makearray(2,xdim)); nprotect++;
setrownames(newstates,Xnames,2);
ss = REAL(x);
st = REAL(newstates);
// create storage for new parameters
- if (!is_op) {
- xdim[0] = npars; xdim[1] = nreps;
+ if (do_par_resamp) {
+ xdim[0] = npars; xdim[1] = np;
PROTECT(newparams = makearray(2,xdim)); nprotect++;
setrownames(newparams,Pnames,2);
ps = REAL(params);
@@ -237,11 +239,11 @@
}
// resample
- nosort_resamp(nreps,REAL(weights),sample,0);
- for (k = 0; k < nreps; k++) { // copy the particles
+ nosort_resamp(nreps,REAL(weights),np,sample,0);
+ for (k = 0; k < np; k++) { // copy the particles
for (j = 0, xx = ss+nvars*sample[k]; j < nvars; j++, st++, xx++)
*st = *xx;
- if (!is_op) {
+ if (do_par_resamp) {
for (j = 0, xp = ps+npars*sample[k]; j < npars; j++, pt++, xp++)
*pt = *xp;
}
@@ -259,7 +261,9 @@
if (do_rw) { // if random walk, adjust prediction variance and move particles
xx = REAL(rw_sd);
- xp = (all_fail&&(!is_op)) ? REAL(params) : REAL(newparams);
+ xp = (all_fail || !do_par_resamp) ? REAL(params) : REAL(newparams);
+ nreps = (all_fail) ? nreps : np;
+
for (j = 0; j < nrw; j++) {
offset = pidx[j];
vsq = xx[j];
@@ -295,7 +299,7 @@
SET_ELEMENT(retval,3,newstates);
}
- if (all_fail||is_op) {
+ if (all_fail || !do_par_resamp) {
SET_ELEMENT(retval,4,params);
} else {
SET_ELEMENT(retval,4,newparams);
Modified: pkg/tests/ou2-forecast.R
===================================================================
--- pkg/tests/ou2-forecast.R 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/ou2-forecast.R 2011-05-20 12:19:27 UTC (rev 490)
@@ -18,7 +18,7 @@
mse <- array(dim=c(2,9,10))
t0 <- seq(from=10,to=90,by=10)
for (k in seq_along(t0)) {
- pp[c("x1.0","x2.0"),] <- pf$saved.states[c("x1","x2"),,tm==t0[k]]
+ pp[c("x1.0","x2.0"),] <- pf$saved.states[tm==t0[k]][[1]][c("x1","x2"),]
inds <- which(tm>t0[k]&tm<=t0[k]+10)
Y <- simulate(ou2,params=pp,obs=TRUE,t0=t0[k],times=tm[inds])
mn <- apply(Y,c(1,3),mean)
Modified: pkg/tests/ou2-forecast.Rout.save
===================================================================
--- pkg/tests/ou2-forecast.Rout.save 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/ou2-forecast.Rout.save 2011-05-20 12:19:27 UTC (rev 490)
@@ -39,7 +39,7 @@
> mse <- array(dim=c(2,9,10))
> t0 <- seq(from=10,to=90,by=10)
> for (k in seq_along(t0)) {
-+ pp[c("x1.0","x2.0"),] <- pf$saved.states[c("x1","x2"),,tm==t0[k]]
++ pp[c("x1.0","x2.0"),] <- pf$saved.states[tm==t0[k]][[1]][c("x1","x2"),]
+ inds <- which(tm>t0[k]&tm<=t0[k]+10)
+ Y <- simulate(ou2,params=pp,obs=TRUE,t0=t0[k],times=tm[inds])
+ mn <- apply(Y,c(1,3),mean)
Modified: pkg/tests/ou2-mif.R
===================================================================
--- pkg/tests/ou2-mif.R 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/ou2-mif.R 2011-05-20 12:19:27 UTC (rev 490)
@@ -150,7 +150,8 @@
Nmif=2,
ivps=c("x1.0","x2.0"),
rw.sd=c(x1.0=5,x2.0=5,alpha.2=0.1,alpha.3=0.2),
- Np=1000,cooling.factor=0.95,ic.lag=10,var.factor=1
+ Np=function(k)if(k<10) 2000 else 500,
+ cooling.factor=0.95,ic.lag=10,var.factor=1
)
fit <- continue(fit)
fit <- continue(fit,Nmif=2)
Modified: pkg/tests/ou2-mif.Rout.save
===================================================================
--- pkg/tests/ou2-mif.Rout.save 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/ou2-mif.Rout.save 2011-05-20 12:19:27 UTC (rev 490)
@@ -1,7 +1,8 @@
-R version 2.11.1 (2010-05-31)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.12.2 (2011-02-25)
+Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
+Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
@@ -16,6 +17,9 @@
Type 'q()' to quit R.
> library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
>
> data(ou2)
>
@@ -138,7 +142,8 @@
+ Np=-10,cooling.factor=0.95,ic.lag=10,var.factor=1
+ )
+ )
-Error : mif error: 'Np' must be a positive integer
+Error in mif.internal(object = object, Nmif = Nmif, start = start, pars = pars, :
+ number of particles, 'Np', must always be positive
>
> try(
+ mif(
@@ -188,7 +193,8 @@
+ Nmif=2,
+ ivps=c("x1.0","x2.0"),
+ rw.sd=c(x1.0=5,x2.0=5,alpha.2=0.1,alpha.3=0.2),
-+ Np=1000,cooling.factor=0.95,ic.lag=10,var.factor=1
++ Np=function(k)if(k<10) 2000 else 500,
++ cooling.factor=0.95,ic.lag=10,var.factor=1
+ )
> fit <- continue(fit)
> fit <- continue(fit,Nmif=2)
Modified: pkg/tests/ou2-procmeas.R
===================================================================
--- pkg/tests/ou2-procmeas.R 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/ou2-procmeas.R 2011-05-20 12:19:27 UTC (rev 490)
@@ -2,13 +2,15 @@
data(ou2)
+po <- window(ou2,end=10)
+
set.seed(3434388L)
-pmat <- do.call(cbind,rep(list(coef(ou2)),3))
-sims <- simulate(ou2,states=T,obs=T,params=pmat)
+pmat <- do.call(cbind,rep(list(coef(po)),3))
+sims <- simulate(po,states=T,obs=T,params=pmat)
-dp <- dprocess(ou2,x=sims$states,times=time(ou2),params=pmat,log=T)
-dm <- dmeasure(ou2,x=sims$states,y=obs(ou2),times=time(ou2),params=pmat,log=T)
+dp <- dprocess(po,x=sims$states,times=time(po),params=pmat,log=T)
+dm <- dmeasure(po,x=sims$states,y=obs(po),times=time(po),params=pmat,log=T)
apply(dp,1,sum)
apply(dm,1,sum)
Modified: pkg/tests/ou2-procmeas.Rout.save
===================================================================
--- pkg/tests/ou2-procmeas.Rout.save 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/ou2-procmeas.Rout.save 2011-05-20 12:19:27 UTC (rev 490)
@@ -23,16 +23,18 @@
>
> data(ou2)
>
+> po <- window(ou2,end=10)
+>
> set.seed(3434388L)
>
-> pmat <- do.call(cbind,rep(list(coef(ou2)),3))
-> sims <- simulate(ou2,states=T,obs=T,params=pmat)
+> pmat <- do.call(cbind,rep(list(coef(po)),3))
+> sims <- simulate(po,states=T,obs=T,params=pmat)
>
-> dp <- dprocess(ou2,x=sims$states,times=time(ou2),params=pmat,log=T)
-> dm <- dmeasure(ou2,x=sims$states,y=obs(ou2),times=time(ou2),params=pmat,log=T)
+> dp <- dprocess(po,x=sims$states,times=time(po),params=pmat,log=T)
+> dm <- dmeasure(po,x=sims$states,y=obs(po),times=time(po),params=pmat,log=T)
>
> apply(dp,1,sum)
-[1] -469.8877 -451.1436 -445.4460
+[1] -47.20607 -37.35755 -46.26917
> apply(dm,1,sum)
-[1] -16202.19 -12225.52 -12082.81
+[1] -850.6711 -616.1925 -993.1957
>
Modified: pkg/tests/ou2-trajmatch.Rout.save
===================================================================
--- pkg/tests/ou2-trajmatch.Rout.save 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/ou2-trajmatch.Rout.save 2011-05-20 12:19:27 UTC (rev 490)
@@ -1,7 +1,8 @@
-R version 2.11.1 (2010-05-31)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.12.2 (2011-02-25)
+Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
+Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
@@ -43,7 +44,7 @@
> range(x['conv',])
[1] 0 0
> range(x['loglik',])
-[1] -688.0077 -567.8361
+[1] -687.0136 -614.8054
> print(
+ cbind(
+ truth=true.p,
@@ -54,16 +55,16 @@
+ digits=4
+ )
truth mean.est bias se
-alpha.1 0.8 0.6724 -0.12758 0.3949
+alpha.1 0.8 0.7082 -0.09185 0.3813
alpha.2 -0.5 -0.5000 0.00000 0.0000
alpha.3 0.3 0.3000 0.00000 0.0000
-alpha.4 0.9 0.8318 -0.06822 0.4922
+alpha.4 0.9 0.9704 0.07039 0.2009
sigma.1 3.0 3.0000 0.00000 0.0000
sigma.2 -0.5 -0.5000 0.00000 0.0000
sigma.3 2.0 2.0000 0.00000 0.0000
-tau 1.0 6.0282 5.02824 1.4802
-x1.0 -3.0 -3.3701 -0.37012 2.7226
-x2.0 4.0 2.3504 -1.64964 4.4465
+tau 1.0 6.0401 5.04009 0.8810
+x1.0 -3.0 -4.2514 -1.25144 2.0408
+x2.0 4.0 4.2243 0.22427 6.3684
>
> summary(traj.match(ou2,est=c('alpha.1','alpha.4','x1.0','x2.0','tau'),method="subplex",maxit=100))
$params
Modified: pkg/tests/pfilter.R
===================================================================
--- pkg/tests/pfilter.R 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/pfilter.R 2011-05-20 12:19:27 UTC (rev 490)
@@ -2,11 +2,20 @@
data(ou2)
+set.seed(9994847L)
+
pf <- pfilter(ou2,Np=1000,seed=343439L)
print(coef(ou2,c('x1.0','x2.0','alpha.1','alpha.4')),digits=4)
cat("particle filter log likelihood at truth\n")
print(pf$loglik,digits=4)
+pf <- replicate(n=10,pfilter(ou2,Np=function(k)if(k<10) 10000 else 500))
+pf.ll <- sapply(pf,logLik)
+ll.est <- log(mean(exp(pf.ll-mean(pf.ll))))+mean(pf.ll)
+ll.se <- sd(exp(pf.ll-mean(pf.ll)))/exp(ll.est-mean(pf.ll))/sqrt(length(pf))
+print(round(c(loglik=ll.est,loglik.se=ll.se),digits=2))
+
data(euler.sir)
pf <- pfilter(euler.sir,Np=200,seed=394343L)
print(pf$loglik,digits=4)
+
Modified: pkg/tests/pfilter.Rout.save
===================================================================
--- pkg/tests/pfilter.Rout.save 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/pfilter.Rout.save 2011-05-20 12:19:27 UTC (rev 490)
@@ -1,7 +1,8 @@
-R version 2.11.1 (2010-05-31)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.12.2 (2011-02-25)
+Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
+Platform: x86_64-unknown-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
@@ -16,9 +17,14 @@
Type 'q()' to quit R.
> library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
>
> data(ou2)
>
+> set.seed(9994847L)
+>
> pf <- pfilter(ou2,Np=1000,seed=343439L)
> print(coef(ou2,c('x1.0','x2.0','alpha.1','alpha.4')),digits=4)
x1.0 x2.0 alpha.1 alpha.4
@@ -28,8 +34,17 @@
> print(pf$loglik,digits=4)
[1] -476.6
>
+> pf <- replicate(n=10,pfilter(ou2,Np=function(k)if(k<10) 10000 else 500))
+> pf.ll <- sapply(pf,logLik)
+> ll.est <- log(mean(exp(pf.ll-mean(pf.ll))))+mean(pf.ll)
+> ll.se <- sd(exp(pf.ll-mean(pf.ll)))/exp(ll.est-mean(pf.ll))/sqrt(length(pf))
+> print(round(c(loglik=ll.est,loglik.se=ll.se),digits=2))
+ loglik loglik.se
+ -479.61 0.46
+>
> data(euler.sir)
> pf <- pfilter(euler.sir,Np=200,seed=394343L)
> print(pf$loglik,digits=4)
-[1] -872.3
+[1] -884.2
>
+>
Modified: pkg/tests/sir-icfit.R
===================================================================
--- pkg/tests/sir-icfit.R 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/sir-icfit.R 2011-05-20 12:19:27 UTC (rev 490)
@@ -116,7 +116,7 @@
}
for (k in seq_len(5)) {
pf <- pfilter(po,params=pp,save.params=TRUE)
- pp <- pf at saved.params[,,nd]
+ pp <- pf at saved.params[[nd]]
}
guesses <- sobol.design(
lower=apply(pp[est,],1,min),
Modified: pkg/tests/sir-icfit.Rout.save
===================================================================
--- pkg/tests/sir-icfit.Rout.save 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/sir-icfit.Rout.save 2011-05-20 12:19:27 UTC (rev 490)
@@ -669,7 +669,7 @@
+ }
+ for (k in seq_len(5)) {
+ pf <- pfilter(po,params=pp,save.params=TRUE)
-+ pp <- pf at saved.params[,,nd]
++ pp <- pf at saved.params[[nd]]
+ }
+ guesses <- sobol.design(
+ lower=apply(pp[est,],1,min),
Modified: pkg/tests/synlik.R
===================================================================
--- pkg/tests/synlik.R 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/synlik.R 2011-05-20 12:19:27 UTC (rev 490)
@@ -4,12 +4,71 @@
set.seed(6457673L)
+kalman.filter <- function (y, x0, a, b, sigma, tau) {
+ n <- nrow(y)
+ ntimes <- ncol(y)
+ sigma.sq <- sigma%*%t(sigma)
+ tau.sq <- tau%*%t(tau)
+ inv.tau.sq <- solve(tau.sq)
+ cond.dev <- numeric(ntimes)
+ filter.mean <- matrix(0,n,ntimes)
+ pred.mean <- matrix(0,n,ntimes)
+ pred.var <- array(0,dim=c(n,n,ntimes))
+ dev <- 0
+ m <- x0
+ v <- diag(0,n)
+ for (k in seq(length=ntimes)) {
+ pred.mean[,k] <- M <- a%*%m
+ pred.var[,,k] <- V <- a%*%v%*%t(a)+sigma.sq
+ q <- b%*%V%*%t(b)+tau.sq
+ r <- y[,k]-b%*%M
+ cond.dev[k] <- n*log(2*pi)+log(det(q))+t(r)%*%solve(q,r)
+ dev <- dev+cond.dev[k]
+ q <- t(b)%*%inv.tau.sq%*%b+solve(V)
+ v <- solve(q)
+ filter.mean[,k] <- m <- v%*%(t(b)%*%inv.tau.sq%*%y[,k]+solve(V,M))
+ }
+ list(
+ pred.mean=pred.mean,
+ pred.var=pred.var,
+ filter.mean=filter.mean,
+ cond.loglik=-0.5*cond.dev,
+ loglik=-0.5*dev
+ )
+}
+
+kalman <- function (x, object, params) {
+ y <- data.array(object)
+ p <- params
+ p[names(x)] <- x
+ x0 <- init.state(object,params=p)
+ a <- matrix(p[c('alpha.1','alpha.2','alpha.3','alpha.4')],2,2)
+ b <- diag(1,2)
+ sigma <- matrix(p[c('sigma.1','sigma.2','sigma.2','sigma.3')],2,2)
+ sigma[1,2] <- 0
+ tau <- diag(p['tau'],2,2)
+ -kalman.filter(y,x0,a,b,sigma,tau)$loglik
+}
+
po <- window(ou2,end=5)
-pb <- replicate(n=100,logLik(probe(po,nsim=1000,probes=function(x)x)))
-pf <- replicate(n=100,logLik(pfilter(po,Np=1000)))
+# exact likelihood
+p.truth <- coef(po)
+loglik.truth <- -kalman(p.truth,po,p.truth)
-kruskal.test(list(pb,pf))
-ks.test(pf,pb)
-qqplot(pf,pb)
+## likelihood from probes (works since ou2 is Gaussian)
+loglik.probe <- replicate(n=500,logLik(probe(po,nsim=200,probes=function(x)x)))
+## likelihood from particle filters
+loglik.pfilter <- replicate(n=500,logLik(pfilter(po,Np=200)))
+
+kruskal.test(list(loglik.probe,loglik.pfilter))
+wilcox.test(loglik.probe,loglik.pfilter)
+ks.test(loglik.pfilter,loglik.probe)
+qqplot(loglik.pfilter,loglik.probe)
abline(a=0,b=1)
+abline(v=loglik.truth)
+
+##require(ggplot2)
+## x <- data.frame(probe=loglik.probe,pfilte=loglik.pfilter)
+## ggplot(data=melt(x))+geom_bar(aes(x=value,fill=variable),position="dodge")
+
Modified: pkg/tests/synlik.Rout.save
===================================================================
--- pkg/tests/synlik.Rout.save 2011-05-19 21:51:10 UTC (rev 489)
+++ pkg/tests/synlik.Rout.save 2011-05-20 12:19:27 UTC (rev 490)
@@ -25,26 +25,92 @@
>
> set.seed(6457673L)
>
+> kalman.filter <- function (y, x0, a, b, sigma, tau) {
++ n <- nrow(y)
++ ntimes <- ncol(y)
++ sigma.sq <- sigma%*%t(sigma)
++ tau.sq <- tau%*%t(tau)
++ inv.tau.sq <- solve(tau.sq)
++ cond.dev <- numeric(ntimes)
++ filter.mean <- matrix(0,n,ntimes)
++ pred.mean <- matrix(0,n,ntimes)
++ pred.var <- array(0,dim=c(n,n,ntimes))
++ dev <- 0
++ m <- x0
++ v <- diag(0,n)
++ for (k in seq(length=ntimes)) {
++ pred.mean[,k] <- M <- a%*%m
++ pred.var[,,k] <- V <- a%*%v%*%t(a)+sigma.sq
++ q <- b%*%V%*%t(b)+tau.sq
++ r <- y[,k]-b%*%M
++ cond.dev[k] <- n*log(2*pi)+log(det(q))+t(r)%*%solve(q,r)
++ dev <- dev+cond.dev[k]
++ q <- t(b)%*%inv.tau.sq%*%b+solve(V)
++ v <- solve(q)
++ filter.mean[,k] <- m <- v%*%(t(b)%*%inv.tau.sq%*%y[,k]+solve(V,M))
++ }
++ list(
++ pred.mean=pred.mean,
++ pred.var=pred.var,
++ filter.mean=filter.mean,
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 490
More information about the pomp-commits
mailing list