[Pomp-commits] r315 - in pkg: . R inst inst/doc src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Sep 20 17:48:20 CEST 2010
Author: kingaa
Date: 2010-09-20 17:48:19 +0200 (Mon, 20 Sep 2010)
New Revision: 315
Added:
pkg/src/pfilter.c
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/pfilter.R
pkg/inst/ChangeLog
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/tests/ricker.Rout.save
Log:
- move most time-consuming computations in pfilter into C
- thus: new file src/pfilter.c
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-08-31 20:40:14 UTC (rev 314)
+++ pkg/DESCRIPTION 2010-09-20 15:48:19 UTC (rev 315)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.32-5
-Date: 2010-08-30
+Version: 0.32-6
+Date: 2010-09-20
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>
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-08-31 20:40:14 UTC (rev 314)
+++ pkg/NAMESPACE 2010-09-20 15:48:19 UTC (rev 315)
@@ -9,6 +9,7 @@
SSA_simulator,
R_Euler_Multinom,
D_Euler_Multinom,
+ pfilter_computations,
do_rprocess,
do_dprocess,
do_rmeasure,
Modified: pkg/R/pfilter.R
===================================================================
--- pkg/R/pfilter.R 2010-08-31 20:40:14 UTC (rev 314)
+++ pkg/R/pfilter.R 2010-09-20 15:48:19 UTC (rev 315)
@@ -77,44 +77,48 @@
}
loglik <- rep(NA,ntimes)
- eff.sample.size <- rep(NA,ntimes)
+ eff.sample.size <- numeric(ntimes)
nfail <- 0
npars <- length(rw.names)
+
+ pred.m <- NULL
+ pred.v <- NULL
+ filt.m <- NULL
+
+ ## set up storage for prediction means, variances, etc.
+ if (pred.mean)
+ pred.m <- matrix(
+ data=0,
+ nrow=nvars+npars,
+ ncol=ntimes,
+ dimnames=list(c(statenames,rw.names),NULL)
+ )
- pred.m <- if (pred.mean)
- matrix(
- data=0,
- nrow=nvars+npars,
- ncol=ntimes,
- dimnames=list(c(statenames,rw.names),NULL)
- )
- else NULL
+ if (pred.var)
+ pred.v <- matrix(
+ data=0,
+ nrow=nvars+npars,
+ ncol=ntimes,
+ dimnames=list(c(statenames,rw.names),NULL)
+ )
- pred.v <- if (pred.var)
- matrix(
- data=0,
- nrow=nvars+npars,
- ncol=ntimes,
- dimnames=list(c(statenames,rw.names),NULL)
- )
- else NULL
-
- filt.m <- if (filter.mean)
- if (random.walk)
- matrix(
- data=0,
- nrow=nvars+length(paramnames),
- ncol=ntimes,
- dimnames=list(c(statenames,paramnames),NULL)
- )
- else
- matrix(
- data=0,
- nrow=nvars,
- ncol=ntimes,
- dimnames=list(statenames,NULL)
- )
- else NULL
+ if (filter.mean) {
+ if (random.walk) {
+ filt.m <- matrix(
+ data=0,
+ nrow=nvars+length(paramnames),
+ ncol=ntimes,
+ dimnames=list(c(statenames,paramnames),NULL)
+ )
+ } else {
+ filt.m <- matrix(
+ data=0,
+ nrow=nvars,
+ ncol=ntimes,
+ dimnames=list(statenames,NULL)
+ )
+ }
+ }
for (nt in seq_len(ntimes)) {
@@ -133,33 +137,16 @@
x[,] <- X # ditch the third dimension
- ## prediction means
- if (pred.mean) {
- xx <- try(
- c(
- rowMeans(x),
- rowMeans(params[rw.names,,drop=FALSE])
- ),
- silent=FALSE
- )
- if (inherits(xx,'try-error')) {
- stop(sQuote("pfilter")," error: error in prediction mean computation",call.=FALSE)
- } else {
- pred.m[,nt] <- xx
- }
- }
-
- ## prediction variances
- if (pred.var) {
+ if (pred.var) { ## check for nonfinite state variables and parameters
problem.indices <- unique(which(!is.finite(x),arr.ind=TRUE)[,1])
- if (length(problem.indices)>0) {
+ if (length(problem.indices)>0) { # state variables
stop(
sQuote("pfilter")," error: non-finite state variable(s): ",
paste(rownames(x)[problem.indices],collapse=', '),
call.=FALSE
)
}
- if (random.walk) {
+ if (random.walk) { # parameters (need to be checked only if 'random.walk=TRUE')
problem.indices <- unique(which(!is.finite(params[rw.names,,drop=FALSE]),arr.ind=TRUE)[,1])
if (length(problem.indices)>0) {
stop(
@@ -169,20 +156,8 @@
)
}
}
- xx <- try(
- c(
- apply(x,1,var),
- apply(params[rw.names,,drop=FALSE],1,var)
- ),
- silent=FALSE
- )
- if (inherits(xx,'try-error')) {
- stop(sQuote("pfilter")," error: error in prediction variance computation",call.=FALSE)
- } else {
- pred.v[,nt] <- xx
- }
}
-
+
## determine the weights
weights <- try(
dmeasure(
@@ -190,48 +165,49 @@
y=object at data[,nt,drop=FALSE],
x=X,
times=times[nt+1],
- params=params
+ params=params,
+ log=FALSE
),
silent=FALSE
)
if (inherits(weights,'try-error'))
stop(sQuote("pfilter")," error: error in calculation of weights",call.=FALSE)
if (any(is.na(weights))) {
- ## problem.indices <- which(is.na(weights))
stop(sQuote("pfilter")," error: ",sQuote("dmeasure")," returns NA",call.=FALSE)
}
- ## test for failure to filter
- dim(weights) <- NULL
- failures <- weights < tol
- all.fail <- all(failures)
- if (all.fail) { # all particles are lost
+ ## prediction mean, prediction variance, filtering mean, effective sample size, log-likelihood
+ xx <- try(
+ .Call(
+ pfilter_computations,
+ x,params,
+ random.walk,rw.names,
+ pred.mean,pred.var,
+ filter.mean,weights,tol
+ ),
+ silent=FALSE
+ )
+ if (inherits(xx,'try-error')) {
+ stop(sQuote("pfilter")," error: error in prediction mean/variance computation",call.=FALSE)
+ }
+ all.fail <- xx$fail
+ loglik[nt] <- xx$loglik
+ eff.sample.size[nt] <- xx$ess
+
+ if (pred.mean)
+ pred.m[,nt] <- xx$pm
+ if (pred.var)
+ pred.v[,nt] <- xx$pv
+ if (filter.mean)
+ filt.m[,nt] <- xx$fm
+
+ if (all.fail) { ## all particles are lost
+ nfail <- nfail+1
if (verbose)
message("filtering failure at time t = ",times[nt+1])
- nfail <- nfail+1
if (nfail > max.fail)
stop(sQuote("pfilter")," error: too many filtering failures",call.=FALSE)
- loglik[nt] <- log(tol) # worst log-likelihood
- weights <- rep(1/Np,Np)
- eff.sample.size[nt] <- 0
- } else { # not all particles are lost
- ## compute log-likelihood
- loglik[nt] <- log(mean(weights))
- weights[failures] <- 0
- weights <- weights/sum(weights)
- ## compute effective sample-size
- eff.sample.size[nt] <- 1/(weights%*%weights)
- }
-
- ## compute filtering means
- if (filter.mean) {
- filt.m[statenames,nt] <- x %*% weights
- if (random.walk)
- filt.m[paramnames,nt] <- params %*% weights
- }
-
- ## Matrix with samples (columns) from filtering distribution theta.t | Y.t
- if (!all.fail) {
+ } else { ## matrix with samples (columns) from filtering distribution theta.t | Y.t
sample <- .Call(systematic_resampling,weights)
x <- x[,sample,drop=FALSE]
params <- params[,sample,drop=FALSE]
@@ -240,7 +216,7 @@
## random walk for parameters
if (random.walk) {
pred.v[rw.names,nt] <- pred.v[rw.names,nt]+sigma^2
- params[rw.names,] <- params[rw.names,]+rnorm(n=Np*length(sigma),mean=0,sd=sigma)
+ params[rw.names,] <- rnorm(n=Np*length(sigma),mean=params[rw.names,],sd=sigma)
}
if (save.states) {
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-08-31 20:40:14 UTC (rev 314)
+++ pkg/inst/ChangeLog 2010-09-20 15:48:19 UTC (rev 315)
@@ -1,5 +1,11 @@
+2010-08-31 kingaa
+
+ * [r314] man/pomp.Rd: - minor clarification
+
2010-08-30 kingaa
+ * [r313] pomp:
+ * [r312] inst/ChangeLog: - update ChangeLog
* [r311] DESCRIPTION, NAMESPACE, R, data, inst, man,
pomp/DESCRIPTION, pomp/NAMESPACE, pomp/R, pomp/data, pomp/inst,
pomp/man, pomp/src, pomp/tests, src, tests: - put everything in
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Added: pkg/src/pfilter.c
===================================================================
--- pkg/src/pfilter.c (rev 0)
+++ pkg/src/pfilter.c 2010-09-20 15:48:19 UTC (rev 315)
@@ -0,0 +1,185 @@
+// -*- C++ -*-
+
+#include "pomp_internal.h"
+#include <Rdefines.h>
+
+// examines weights for filtering failure
+// computes log likelihood and effective sample size
+// computes (if desired) prediction mean, prediction variance, filtering mean.
+// it is assumed that ncol(x) == ncol(params).
+// 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 rw, SEXP rw_names,
+ SEXP predmean, SEXP predvar,
+ SEXP filtmean, SEXP weights, SEXP tol)
+{
+ int nprotect = 0;
+ SEXP pm, pv, fm, ess, fail, loglik;
+ SEXP retval, retvalnames;
+ double *xpm, *xpv, *xfm, *xw;
+ SEXP dimX, Pnames, pindex;
+ double *xx, *xp;
+ int *dim, *pidx, lv;
+ int nvars, npars, nrw, nreps, offset, nlost;
+ int do_rw, do_pm, do_pv, do_fm, all_fail = 0;
+ double sum, sumsq, vsq, ws, w, toler;
+ int j, k;
+
+ PROTECT(dimX = GET_DIM(x)); nprotect++;
+ dim = INTEGER(dimX);
+ nvars = dim[0]; nreps = dim[1];
+ xx = REAL(x);
+ do_rw = *(LOGICAL(AS_LOGICAL(rw)));
+ do_pm = *(LOGICAL(AS_LOGICAL(predmean)));
+ do_pv = *(LOGICAL(AS_LOGICAL(predvar)));
+ do_fm = *(LOGICAL(AS_LOGICAL(filtmean)));
+
+ PROTECT(ess = NEW_NUMERIC(1)); nprotect++;
+ PROTECT(loglik = NEW_NUMERIC(1)); nprotect++;
+ PROTECT(fail = NEW_LOGICAL(1)); nprotect++;
+
+ xw = REAL(weights);
+ toler = *(REAL(tol));
+
+ // check the weights and compute sum and sum of squares
+ for (k = 0, w = 0, ws = 0, nlost = 0; k < nreps; k++) {
+ if (xw[k] > toler) {
+ w += xw[k];
+ ws += xw[k]*xw[k];
+ } else { // this particle is lost
+ xw[k] = 0;
+ nlost++;
+ }
+ }
+ if (nlost >= nreps) all_fail = 1; // all particles are lost
+ if (all_fail) {
+ *(REAL(loglik)) = log(toler); // minimum log-likelihood
+ *(REAL(ess)) = 0; // zero effective sample size
+ } else {
+ *(REAL(loglik)) = log(w/((double) nreps)); // mean of weights is likelihood
+ *(REAL(ess)) = w*w/ws; // effective sample size
+ }
+ *(LOGICAL(fail)) = all_fail;
+
+ if (do_rw) {
+ PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
+ PROTECT(pindex = matchnames(Pnames,rw_names)); nprotect++;
+ xp = REAL(params);
+ pidx = INTEGER(pindex);
+ npars = LENGTH(Pnames);
+ nrw = LENGTH(rw_names);
+ lv = nvars+nrw;
+ } else {
+ pidx = NULL;
+ lv = nvars;
+ }
+
+ if (do_pm || do_pv) {
+ PROTECT(pm = NEW_NUMERIC(lv)); nprotect++;
+ xpm = REAL(pm);
+ }
+
+ if (do_pv) {
+ PROTECT(pv = NEW_NUMERIC(lv)); nprotect++;
+ xpv = REAL(pv);
+ }
+
+ if (do_fm) {
+ PROTECT(fm = NEW_NUMERIC(nvars+npars)); nprotect++;
+ xfm = REAL(fm);
+ }
+
+ for (j = 0; j < nvars; j++) { // state variables
+
+ // compute prediction mean
+ if (do_pm || do_pv) {
+ for (k = 0, sum = 0; k < nreps; k++) sum += xx[j+k*nvars];
+ sum /= ((double) nreps);
+ xpm[j] = sum;
+ }
+
+ // compute prediction variance
+ if (do_pv) {
+ for (k = 0, sumsq = 0; k < nreps; k++) {
+ vsq = xx[j+k*nvars]-sum;
+ sumsq += vsq*vsq;
+ }
+ xpv[j] = sumsq / ((double) (nreps - 1));
+ }
+
+ // compute filter mean
+ if (do_fm) {
+ if (all_fail) {
+ for (k = 0, ws = 0; k < nreps; k++) ws += xx[j+k*nvars];
+ xfm[j] = ws/((double) nreps);
+ } else {
+ for (k = 0, ws = 0; k < nreps; k++) ws += xx[j+k*nvars]*xw[k];
+ xfm[j] = ws/w;
+ }
+ }
+
+ }
+
+ // compute means and variances for parameters
+ if (do_rw) {
+ for (j = 0; j < nrw; j++) {
+ offset = pidx[j];
+
+ if (do_pm || do_pv) {
+ for (k = 0, sum = 0; k < nreps; k++) sum += xp[offset+k*npars];
+ sum /= ((double) nreps);
+ xpm[nvars+j] = sum;
+ }
+
+ if (do_pv) {
+ for (k = 0, sumsq = 0; k < nreps; k++) {
+ vsq = xp[offset+k*npars]-sum;
+ sumsq += vsq*vsq;
+ }
+ xpv[nvars+j] = sumsq / ((double) (nreps - 1));
+ }
+
+ }
+ if (do_fm) {
+ for (j = 0; j < npars; j++) {
+ if (all_fail) {
+ for (k = 0, ws = 0; k < nreps; k++) ws += xp[j+k*npars];
+ xfm[nvars+j] = ws/((double) nreps);
+ } else {
+ for (k = 0, ws = 0; k < nreps; k++) ws += xp[j+k*npars]*xw[k];
+ xfm[nvars+j] = ws/w;
+ }
+ }
+ }
+ }
+
+ PROTECT(retval = NEW_LIST(6)); nprotect++;
+ PROTECT(retvalnames = NEW_CHARACTER(6)); nprotect++;
+ SET_STRING_ELT(retvalnames,0,mkChar("fail"));
+ SET_STRING_ELT(retvalnames,1,mkChar("loglik"));
+ SET_STRING_ELT(retvalnames,2,mkChar("ess"));
+ SET_STRING_ELT(retvalnames,3,mkChar("pm"));
+ SET_STRING_ELT(retvalnames,4,mkChar("pv"));
+ SET_STRING_ELT(retvalnames,5,mkChar("fm"));
+ SET_NAMES(retval,retvalnames);
+
+ SET_ELEMENT(retval,0,fail);
+ SET_ELEMENT(retval,1,loglik);
+ SET_ELEMENT(retval,2,ess);
+
+ if (do_pm) {
+ SET_ELEMENT(retval,3,pm);
+ }
+ if (do_pv) {
+ SET_ELEMENT(retval,4,pv);
+ }
+ if (do_fm) {
+ SET_ELEMENT(retval,5,fm);
+ }
+
+ UNPROTECT(nprotect);
+ return(retval);
+}
+
Property changes on: pkg/src/pfilter.c
___________________________________________________________________
Added: svn:eol-style
+ native
Modified: pkg/tests/ricker.Rout.save
===================================================================
--- pkg/tests/ricker.Rout.save 2010-08-31 20:40:14 UTC (rev 314)
+++ pkg/tests/ricker.Rout.save 2010-09-20 15:48:19 UTC (rev 315)
@@ -47,12 +47,12 @@
>
> print(res,digits=3)
truth MLE guess
-log.r 3.80 3.73 0.693
-log.sigma -1.20 -1.38 0.000
-log.phi 2.30 2.31 2.996
+log.r 3.80 3.79 0.693
+log.sigma -1.20 -1.40 0.000
+log.phi 2.30 2.29 2.996
N.0 7.00 7.00 7.000
e.0 0.00 0.00 0.000
-loglik -142.74 -143.25 -266.548
+loglik -142.74 -142.24 -266.548
>
> tj.1 <- trajectory(ricker)
> plot(time(ricker),tj.1[1,,-1],type='l')
More information about the pomp-commits
mailing list