[Pomp-commits] r352 - in pkg: . R inst inst/doc src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Sep 30 00:05:47 CEST 2010
Author: kingaa
Date: 2010-09-30 00:05:47 +0200 (Thu, 30 Sep 2010)
New Revision: 352
Added:
pkg/src/probe_ccf.c
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/basic-probes.R
pkg/inst/ChangeLog
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/tests/ou2-probe.R
Log:
- add in a new probe, 'probe.ccf', which computes the CCF of two variables at various lags
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-09-29 21:06:10 UTC (rev 351)
+++ pkg/DESCRIPTION 2010-09-29 22:05:47 UTC (rev 352)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.33-2
-Date: 2010-09-28
+Version: 0.33-3
+Date: 2010-09-29
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-09-29 21:06:10 UTC (rev 351)
+++ pkg/NAMESPACE 2010-09-29 22:05:47 UTC (rev 352)
@@ -15,6 +15,7 @@
probe_marginal_setup,
probe_marginal_solve,
probe_acf,
+ probe_ccf,
probe_nlar,
do_rprocess,
do_dprocess,
@@ -73,6 +74,7 @@
probe.period,
probe.quantile,
probe.acf,
+ probe.ccf,
probe.nlar,
probe.cov,
probe.cor,
Modified: pkg/R/basic-probes.R
===================================================================
--- pkg/R/basic-probes.R 2010-09-29 21:06:10 UTC (rev 351)
+++ pkg/R/basic-probes.R 2010-09-29 22:05:47 UTC (rev 352)
@@ -120,6 +120,18 @@
)
}
+probe.ccf <- function (vars, lags, transform = identity) {
+ transform <- match.fun(transform)
+ if (length(vars)!=2)
+ stop(sQuote("vars")," must name two variables")
+ function (y) .Call(
+ probe_ccf,
+ x=transform(y[vars[1],,drop=TRUE]),
+ y=transform(y[vars[2],,drop=TRUE]),
+ lags=lags
+ )
+}
+
probe.marginal <- function (var, ref, order = 3, diff = 1, transform = identity) {
if (length(var)>1) stop(sQuote("probe.marginal")," is a univariate probe")
transform <- match.fun(transform)
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-09-29 21:06:10 UTC (rev 351)
+++ pkg/inst/ChangeLog 2010-09-29 22:05:47 UTC (rev 352)
@@ -1,5 +1,60 @@
+2010-09-29 kingaa
+
+ * [r351] R/basic-probes.R, src/mat.c, src/pomp_internal.h,
+ src/pomp_mat.h, src/probe_marginal.c, src/probe_nlar.c: - remove
+ 'mat.c' in favor of inline functions for regression, defined in
+ the new header 'pomp_mat.h'
+ * [r350] src/euler.c, src/mat.c, src/pfilter.c, src/probe.c,
+ src/probe_acf.c, src/probe_marginal.c, src/probe_nlar.c: - remove
+ the last usages of dynamic allocation from 'mat.c'
+ - changes in the compiled codes to eliminate pedantic compiler
+ warnings
+ * [r349] tests/ou2-probe.R, tests/ou2-probe.Rout.save: - an exact
+ check for the correctness of probe.acf is added
+ * [r348] man/pomp-methods.Rd: - keep up with changes
+ * [r347] tests/ou2-probe-match.R, tests/ou2-probe-match.Rout.save,
+ tests/ou2-probe.R, tests/ou2-probe.Rout.save: - change name of
+ file
+ - add in a test of probe.acf and probe.nlar
+ * [r346] src/bspline.c: - use R_alloc instead of Calloc/Free
+ * [r345] src/euler.c, src/rmeasure.c: - use R_alloc instead of
+ Calloc/Free
+ * [r344] NAMESPACE, R/pomp-methods.R: - export the new S3
+ as.data.frame method for pomps
+ * [r343] man/basic-probes.Rd, tests/ricker-probe.Rout.save: -
+ update documentation with an entry for 'probe.nlar'
+
+2010-09-28 kingaa
+
+ * [r339] NAMESPACE, R/basic-probes.R, src/probe_nlar.c,
+ tests/ricker-probe.R: - add new 'probe.nlar': this probe performs
+ a polynomial autoregression on a univariate series and returns
+ the coefficients. Recommended by Wood (2010).
+ * [r338] src/probe_marginal.c: - change names returned by
+ 'probe.marginal'
+ * [r337] R/probe.R: - fix up labels when many probes are used
+ * [r334] man/sobol.Rd: - fix documentation error
+ * [r333] DESCRIPTION, NAMESPACE, R/profile-design.R,
+ R/slice-design.R, R/slice.R, R/sobol.R, inst/NEWS,
+ man/profile-design.Rd, man/slice-design.Rd, man/slice.Rd,
+ man/sobol.Rd: - rearrange the slice and profile design functions
+ to make them easier to use
+ * [r330] src/pomp_internal.h: - use R_STATIC macro
+
+2010-09-26 kingaa
+
+ * [r329] tests/ricker-probe.Rout.save: - improve tests for probes
+ * [r328] tests/ricker-probe.R: - improve tests for probes
+ * [r327] src/probe_marginal.c: - bug fix
+
+2010-09-25 kingaa
+
+ * [r326] DESCRIPTION, R/probe-match.R, R/probe.R, src/probe.c,
+ tests/ricker-probe.R: - fix bug with unnamed probes
+
2010-09-24 kingaa
+ * [r325] inst/ChangeLog: - update ChangeLog
* [r324] DESCRIPTION, NAMESPACE, R/basic-probes.R, inst/ChangeLog,
inst/doc/advanced_topics_in_pomp.pdf, inst/doc/intro_to_pomp.pdf,
inst/doc/pomp.bib, man/basic-probes.Rd, src/Makevars, src/mat.c,
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)
Added: pkg/src/probe_ccf.c
===================================================================
--- pkg/src/probe_ccf.c (rev 0)
+++ pkg/src/probe_ccf.c 2010-09-29 22:05:47 UTC (rev 352)
@@ -0,0 +1,92 @@
+// -*- mode: C++ -*-
+
+#include "pomp_internal.h"
+#include <stdio.h>
+
+static void pomp_ccf_compute (double *ccf, double *x, double *y, int n, int *lags, int nlag);
+
+SEXP probe_ccf (SEXP x, SEXP y, SEXP lags) {
+ int nprotect = 0;
+ SEXP ccf, ccf_names, cccf;
+ SEXP X, Y;
+ int nlag, n;
+ int j, k, l;
+ double *p, *p1;
+ char tmp[BUFSIZ], *nm;
+
+ nlag = LENGTH(lags);
+ PROTECT(lags = AS_INTEGER(lags)); nprotect++;
+
+ n = LENGTH(x); // n = # of observations
+ if (n != LENGTH(y))
+ error("'x' and 'y' must have equal lengths");
+
+ PROTECT(X = duplicate(AS_NUMERIC(x))); nprotect++;
+ PROTECT(Y = duplicate(AS_NUMERIC(y))); nprotect++;
+
+ PROTECT(ccf = NEW_NUMERIC(nlag)); nprotect++;
+
+ pomp_ccf_compute(REAL(ccf),REAL(X),REAL(Y),n,INTEGER(lags),nlag);
+
+ PROTECT(ccf_names = NEW_STRING(nlag)); nprotect++;
+ for (k = 0; k < nlag; k++) {
+ snprintf(tmp,BUFSIZ,"ccf.%d",INTEGER(lags)[k]);
+ SET_STRING_ELT(ccf_names,k,mkChar(tmp));
+ }
+ SET_NAMES(ccf,ccf_names);
+
+ UNPROTECT(nprotect);
+ return ccf;
+}
+
+// vectorized routine for CCF calculation
+// we first center the series and then compute means of products
+static void pomp_ccf_compute (double *ccf, double *x, double *y, int n, int *lags, int nlag) {
+ double xx, *p, *p1, *p2;
+ int j, k, lag, ct;
+
+ // first center x
+ for (k = 0, xx = 0, ct = 0, p = x; k < n; k++, p++) {
+ if (R_FINITE(*p)) {
+ xx += *p;
+ ct++;
+ }
+ }
+ if (ct < 1) error("series 'x' has no data");
+ xx /= ct; // mean of x[j]
+ for (k = 0, p = x; k < n; k++, p++)
+ if (R_FINITE(*p)) *p -= xx;
+
+ // now center y
+ for (k = 0, xx = 0, ct = 0, p = y; k < n; k++, p++) {
+ if (R_FINITE(*p)) {
+ xx += *p;
+ ct++;
+ }
+ }
+ if (ct < 1) error("series 'y' has no data");
+ xx /= ct; // mean of y[j]
+ for (k = 0, p = y; k < n; k++, p++)
+ if (R_FINITE(*p)) *p -= xx;
+
+ // compute covariances
+ for (j = 0, p = ccf; j < nlag; j++, p++) { // loop over lags
+ lag = lags[j];
+ if (lag < 0) {
+ for (k = 0, xx = 0, ct = 0, p1 = x-lag, p2 = y; k < n+lag; k++, p1++, p2++)
+ if (R_FINITE(*p1) && R_FINITE(*p1)) {
+ xx += (*p1)*(*p2);
+ ct++;
+ }
+ *p = (ct > 0) ? xx/ct : R_NaReal;
+ } else {
+ for (k = 0, xx = 0, ct = 0, p1 = x, p2 = y+lag; k < n-lag; k++, p1++, p2++)
+ if (R_FINITE(*p1) && R_FINITE(*p1)) {
+ xx += (*p1)*(*p2);
+ ct++;
+ }
+ *p = (ct > 0) ? xx/ct : R_NaReal;
+ }
+ }
+
+}
Property changes on: pkg/src/probe_ccf.c
___________________________________________________________________
Added: svn:eol-style
+ native
Modified: pkg/tests/ou2-probe.R
===================================================================
--- pkg/tests/ou2-probe.R 2010-09-29 21:06:10 UTC (rev 351)
+++ pkg/tests/ou2-probe.R 2010-09-29 22:05:47 UTC (rev 352)
@@ -77,3 +77,28 @@
y2 <- head(x$y2,-1)
z2 <- tail(x$y2,-1)
max(abs(pb at datvals-c(mean(y1*z1)/mean(x$y1^2),mean(y2*z2)/mean(x$y2^2),mean(y1*z1)/mean(y1*y1),mean(y2*z2)/mean(y2*y2))))
+
+po <- simulate(ou2)
+pb <- probe(
+ po,
+ probes=list(
+ probe.acf(var=c("y1"),lag.max=2,type="cov"),
+ probe.ccf(vars=c("y1","y1"),lags=c(0,1,2))
+ ),
+ nsim=1000,
+ seed=1066L
+ )
+plot(pb)
+summary(pb)
+
+pb <- probe(
+ po,
+ probes=probe.ccf(vars=c("y1","y2"),lags=c(-5,-3,1,4,8)),
+ nsim=1000,
+ seed=1066L
+ )
+plot(pb)
+summary(pb)
+
+
+
More information about the pomp-commits
mailing list