[Pomp-commits] r446 - in pkg: . R inst src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Apr 24 15:35:36 CEST 2011
Author: kingaa
Date: 2011-04-24 15:35:35 +0200 (Sun, 24 Apr 2011)
New Revision: 446
Modified:
pkg/DESCRIPTION
pkg/R/pfilter.R
pkg/inst/ChangeLog
pkg/src/resample.c
Log:
- fix bug in 'pfilter' introduced in version 0.32-6 (rev 315): upon filtering failure, state variables were left un-updated. Thanks to Karim Khader for noticing circumstances that led to the notice of this bug.
- minor rearrangement of the codes in resample.c
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2011-04-20 21:04:07 UTC (rev 445)
+++ pkg/DESCRIPTION 2011-04-24 13:35:35 UTC (rev 446)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.36-4
-Date: 2011-04-20
+Version: 0.36-5
+Date: 2011-04-24
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/R/pfilter.R
===================================================================
--- pkg/R/pfilter.R 2011-04-20 21:04:07 UTC (rev 445)
+++ pkg/R/pfilter.R 2011-04-24 13:35:35 UTC (rev 446)
@@ -217,8 +217,9 @@
nfail <- nfail+1
if (verbose)
message("filtering failure at time t = ",times[nt+1])
- if (nfail > max.fail)
+ if (nfail>max.fail)
stop(sQuote("pfilter")," error: too many filtering failures",call.=FALSE)
+ x[,] <- X[,,1,drop=FALSE]
} else { ## matrix with samples (columns) from filtering distribution theta.t | Y.t
sample <- .Call(systematic_resampling,weights)
x[,] <- X[,sample,1,drop=FALSE]
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2011-04-20 21:04:07 UTC (rev 445)
+++ pkg/inst/ChangeLog 2011-04-24 13:35:35 UTC (rev 446)
@@ -1,3 +1,26 @@
+2011-04-20 kingaa
+
+ * [r445] DESCRIPTION, data/dacca.rda, data/euler.sir.rda,
+ data/gillespie.sir.rda, data/gompertz.rda, data/ou2.rda,
+ data/ricker.rda, data/rw2.rda, data/verhulst.rda,
+ inst/data-R/gompertz.R, 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/ou2-multi-mif.rda,
+ inst/doc/ou2-trajmatch.rda, src/gompertz.c: - improvements to the
+ "intro" vignette
+ - fix problem with gompertz measurement model
+ * [r444] DESCRIPTION, R/pomp-methods.R, data/dacca.rda,
+ data/euler.sir.rda, data/gillespie.sir.rda, data/gompertz.rda,
+ data/ou2.rda, data/ricker.rda, data/rw2.rda, data/verhulst.rda,
+ inst/ChangeLog, inst/data-R/gompertz.R,
+ inst/doc/advanced_topics_in_pomp.pdf,
+ 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/pomp.bib, man/gompertz.Rd, src/gompertz.c: - fix bug in
+ 'states' method when 'states' slot is empty
+ - add Gompertz model
+ - edit 'intro_to_pomp' vignette to focus on Gompertz model
+
2011-04-19 kingaa
* [r443] tests/synlik.Rout.save: - add tests/synlik.Rout.save file
Modified: pkg/src/resample.c
===================================================================
--- pkg/src/resample.c 2011-04-20 21:04:07 UTC (rev 445)
+++ pkg/src/resample.c 2011-04-24 13:35:35 UTC (rev 446)
@@ -3,34 +3,91 @@
#include "pomp_internal.h"
#include <Rdefines.h>
-SEXP systematic_resampling (SEXP weights)
+static void nosort_resamp (int n, double *w, int *p, int offset)
{
- int nprotect = 0;
- int n = length(weights);
- double u, du, *w;
- int i, j, *p;
- SEXP perm;
+ int i, j;
+ double du, u;
- PROTECT(perm = NEW_INTEGER(n)); nprotect++;
-
- p = INTEGER(perm);
- w = REAL(weights);
for (j = 1; j < n; j++) w[j] += w[j-1];
- for (j = 0; j < n; j++) w[j] /= w[n-1];
+ if (w[n-1] <= 0.0)
+ error("non-positive sum of weights");
+
GetRNGstate();
- du = 1.0 / ((double) n);
+ du = w[n-1] / ((double) n);
u = runif(-du,0);
PutRNGstate();
for (i = 0, j = 0; j < n; j++) {
u += du;
while (u > w[i]) i++;
- p[j] = i + 1; // must use 1-based indexing for compatibility with R!
+ p[j] = i;
}
+ if (offset) // add offset if needed
+ for (j = 0; j < n; j++) p[j] += offset;
+}
+
+SEXP systematic_resampling (SEXP weights)
+{
+ int nprotect = 0;
+ double u, du, *w;
+ int i, j, *p, n;
+ SEXP perm;
+
+ n = LENGTH(weights);
+ PROTECT(perm = NEW_INTEGER(n)); nprotect++;
+ nosort_resamp(n,REAL(weights),INTEGER(perm),1);
UNPROTECT(nprotect);
-
return(perm);
}
+SEXP pfilter_resample (SEXP weights, SEXP states, SEXP params)
+{
+ int nprotect = 0;
+ int xdim[2], nvars, npars, nreps;
+ SEXP sdim, pdim, newstates, newparams, sample;
+ SEXP retval, retvalnames;
+ int i, j, n;
+ double *ss, *st, *ps, *pt, *xx;
+
+ PROTECT(sdim = GET_DIM(states)); nprotect++;
+ nvars = INTEGER(sdim)[0]; nreps = INTEGER(sdim)[1];
+
+ xdim[0] = nvars; xdim[1] = nreps;
+ PROTECT(newstates = makearray(2,xdim)); nprotect++;
+ setrownames(newstates,GET_ROWNAMES(GET_DIMNAMES(states)),2);
+ ss = REAL(states);
+ st = REAL(newstates);
+
+ PROTECT(pdim = GET_DIM(params)); nprotect++;
+ npars = INTEGER(pdim)[0];
+ if (nreps != INTEGER(pdim)[1])
+ error("'states' and 'params' do not agree in second dimension");
+
+ xdim[0] = npars; xdim[1] = nreps;
+ PROTECT(newparams = makearray(2,xdim)); nprotect++;
+ setrownames(newparams,GET_ROWNAMES(GET_DIMNAMES(params)),2);
+ ps = REAL(params);
+ pt = REAL(newparams);
+
+ {
+ int sample[nreps];
+ nosort_resamp(nreps,REAL(weights),sample,0);
+ for (i = 0; i < nreps; i++) {
+ for (j = 0, xx = ss+nvars*sample[i]; j < nvars; j++, st++, xx++) *st = *xx;
+ for (j = 0, xx = ps+npars*sample[i]; j < npars; j++, pt++, xx++) *pt = *xx;
+ }
+ }
+
+ PROTECT(retval = NEW_LIST(2)); nprotect++;
+ SET_ELEMENT(retval,0,newstates);
+ SET_ELEMENT(retval,1,newparams);
+ PROTECT(retvalnames = NEW_CHARACTER(2)); nprotect++;
+ SET_STRING_ELT(retvalnames,0,mkChar("states"));
+ SET_STRING_ELT(retvalnames,1,mkChar("params"));
+ SET_NAMES(retval,retvalnames);
+
+ UNPROTECT(nprotect);
+ return retval;
+}
More information about the pomp-commits
mailing list