[Pomp-commits] r339 - in pkg: . R src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 28 19:34:28 CEST 2010
Author: kingaa
Date: 2010-09-28 19:34:27 +0200 (Tue, 28 Sep 2010)
New Revision: 339
Added:
pkg/src/probe_nlar.c
Modified:
pkg/NAMESPACE
pkg/R/basic-probes.R
pkg/tests/ricker-probe.R
Log:
- add new 'probe.nlar': this probe performs a polynomial autoregression on a univariate series and returns the coefficients. Recommended by Wood (2010).
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-09-28 17:33:13 UTC (rev 338)
+++ pkg/NAMESPACE 2010-09-28 17:34:27 UTC (rev 339)
@@ -15,6 +15,7 @@
probe_marginal_setup,
probe_marginal_solve,
probe_acf,
+ probe_nlar,
do_rprocess,
do_dprocess,
do_rmeasure,
@@ -71,6 +72,7 @@
probe.period,
probe.quantile,
probe.acf,
+ probe.nlar,
probe.cov,
probe.cor,
probe.marginal,
Modified: pkg/R/basic-probes.R
===================================================================
--- pkg/R/basic-probes.R 2010-09-28 17:33:13 UTC (rev 338)
+++ pkg/R/basic-probes.R 2010-09-28 17:34:27 UTC (rev 339)
@@ -131,3 +131,24 @@
corr=corr
)
}
+
+probe.nlar <- function (var, lags, powers, transform = identity) {
+ transform <- match.fun(transform)
+ if (any(lags<1)||any(powers<1))
+ stop(sQuote("lags")," and ",sQuote("powers")," must be positive integers")
+ if (length(lags)<length(powers)) {
+ if (length(lags)>1) stop(sQuote("lags")," must match ",sQuote("powers")," in length, or have length 1")
+ lags <- rep(lags,length(powers))
+ } else if (length(lags)>length(powers)) {
+ if (length(powers)>1) stop(sQuote("powers")," must match ",sQuote("lags")," in length, or have length 1")
+ powers <- rep(powers,length(lags))
+ }
+ lags <- as.integer(lags)
+ powers <- as.integer(powers)
+ function (y) .Call(
+ probe_nlar,
+ x=transform(y[var,,drop=FALSE]),
+ lags=lags,
+ powers=powers
+ )
+}
Added: pkg/src/probe_nlar.c
===================================================================
--- pkg/src/probe_nlar.c (rev 0)
+++ pkg/src/probe_nlar.c 2010-09-28 17:34:27 UTC (rev 339)
@@ -0,0 +1,139 @@
+// -*- mode: C++; -*-
+
+#include "pomp_internal.h"
+#include <stdio.h>
+
+static void pomp_nlar(double *beta, double *y, int n, int nterms, int *lag, int *power);
+
+SEXP probe_nlar (SEXP x, SEXP lags, SEXP powers) {
+ int nprotect = 0;
+ SEXP y, beta, beta_names;
+ int n, nterms;
+ int k;
+ char tmp[BUFSIZ];
+
+ n = LENGTH(x); // n = # of observations
+ nterms = LENGTH(lags);
+
+ PROTECT(y = duplicate(AS_NUMERIC(x))); nprotect++;
+ PROTECT(beta = NEW_NUMERIC(nterms)); nprotect++;
+
+ pomp_nlar(REAL(beta),REAL(y),n,nterms,INTEGER(lags),INTEGER(powers));
+
+ PROTECT(beta_names = NEW_STRING(nterms)); nprotect++;
+ for (k = 0; k < nterms; k++) {
+ snprintf(tmp,BUFSIZ,"nlar.%ld^%ld",INTEGER(lags)[k],INTEGER(powers)[k]);
+ SET_STRING_ELT(beta_names,k,mkChar(tmp));
+ }
+ SET_NAMES(beta,beta_names);
+
+ UNPROTECT(nprotect);
+ return beta;
+}
+
+// Code for polynomial auto-regression
+// The original version of the following code is due to Simon N. Wood.
+// Modifications by AAK.
+static void pomp_nlar(double *beta, double *y, int n,
+ int nterms, int *lag, int *power) {
+ // 'x' is an n vector of data.
+ // 'nterms' gives the number of terms on the rhs of the autoregression.
+ // 'lag[i]' gives the lag of the ith term on the rhs.
+ // 'power[i]' gives the power to which the ith rhs term should be raised.
+ // 'beta' contains the ar coefficients
+
+ int maxlag, nx, ny, ok, ct;
+ int one = 1;
+ int i, j, k;
+ double xx, *yp;
+ double obs1;
+
+ // find maxlag
+ for (maxlag = 0, i = 0; i < nterms; i++)
+ maxlag = (lag[i] > maxlag) ? lag[i] : maxlag;
+ ny = n - maxlag; // maximum response vector length
+
+ // compute the series mean
+ for (j = 0, ct = 0, xx = 0.0; j < n; j++)
+ if (R_FINITE(y[j])) {
+ xx += y[j];
+ ct++;
+ }
+ xx /= ct; // series mean
+
+ // center the whole series
+ // also check to see if there is any variability in the predictors
+ ok = 0;
+ obs1 = R_NaReal;
+ for (j = 0; j < n; j++) {
+ if (R_FINITE(y[j])) {
+ if (!ok && (j < ny)) { // j < ny means x[j] is a predictor
+ if (!R_FINITE(obs1)) obs1 = y[j];
+ else if (y[j] != obs1) ok = 1;
+ }
+ y[j] -= xx; // subtracting series mean
+ }
+ }
+
+ if (!ok) { // data had no variability
+
+ for (i = 0; i < nterms; i++) beta[i] = 0.0;
+
+ } else { // data not all the same
+
+ double *X, *Xp;
+ int finite[ny];
+
+ // test for NA rows in model matrix and response vector
+ for (nx = 0, yp = y+maxlag, j = 0; j < ny; j++) {
+ finite[j] = (R_FINITE(yp[j])) ? 1 : 0; // finite response?
+ for (i = 0; i < nterms; i++) {
+ finite[j] = (R_FINITE(yp[j-lag[i]])) ? finite[j] : 0; // finite in model matrix row j, column i
+ }
+ if (finite[j]) nx++;
+ }
+ // nx is now the number of non-NA rows in the model matrix
+
+ // build the model matrix, omitting NA rows
+ X = (double *) Calloc(nx*nterms,double);
+ for (Xp = X, i = 0; i < nterms; i++) { // work through the terms
+ for (j = 0; j < ny; j++) {
+ if (finite[j]) {
+ xx = yp[j-lag[i]]; // current predictor
+ *Xp = xx;
+ for (k = 1; k < power[i]; k++) *Xp *= xx; // raise to the appropriate power
+ Xp++;
+ }
+ }
+ }
+ // X is now the nx by nterms model matrix
+
+ // drop the NA rows from the response data
+ for (i = 0, j = 0; j < ny; j++) {
+ if (finite[j]) { // keep this row
+ yp[i++] = yp[j];
+ }
+ }
+ // response vector is now length nx
+
+ {
+ double tau[nterms], b[nterms];
+ int pivot[nterms];
+
+ // QR decompose the model matrix
+ for (i = 0; i < nterms; i++) pivot[i] = 0;
+ pomp_qr(X,&nx,&nterms,pivot,tau);
+ // solve R b = Q'y for b
+ pomp_qrqy(yp,X,tau,&nx,&one,&nterms,&one,&one); // y <- Q'y
+ pomp_backsolve(X,&nx,&nterms,yp,b,&one); // b <- R^{-1} Q'y
+
+ // store b in kth column of beta
+ for (i = 0; i < nterms; i++) beta[pivot[i]] = b[i];
+
+ }
+
+ Free(X);
+
+ }
+
+}
Property changes on: pkg/src/probe_nlar.c
___________________________________________________________________
Added: svn:eol-style
+ native
Modified: pkg/tests/ricker-probe.R
===================================================================
--- pkg/tests/ricker-probe.R 2010-09-28 17:33:13 UTC (rev 338)
+++ pkg/tests/ricker-probe.R 2010-09-28 17:34:27 UTC (rev 339)
@@ -8,6 +8,7 @@
z <- as.numeric(data.array(ricker))
po <- ricker
+
pb <- probe(
po,
probes=probe.median("y"),
@@ -19,6 +20,33 @@
pb <- probe(
po,
+ probes=probe.nlar("y",lags=c(1,2,3),powers=c(1,1,1),transform="sqrt"),
+ nsim=1000,
+ seed=838775L
+ )
+plot(pb)
+summary(pb)
+
+pb <- probe(
+ po,
+ probes=probe.nlar("y",lags=c(1,2,3),powers=1,transform="sqrt"),
+ nsim=1000,
+ seed=838775L
+ )
+plot(pb)
+summary(pb)
+
+pb <- probe(
+ po,
+ probes=probe.nlar("y",lags=1,powers=c(1,2,3),transform="sqrt"),
+ nsim=1000,
+ seed=838775L
+ )
+plot(pb)
+summary(pb)
+
+pb <- probe(
+ po,
probes=probe.marginal(
var="y",
transform=sqrt,
@@ -92,6 +120,37 @@
pb <- probe(
po,
+ probes=probe.nlar(
+ var="y",
+ transform=sqrt,
+ lags=1,
+ powers=c(1,2,3)
+ ),
+ nsim=1000,
+ seed=838775L
+ )
+pb at datvals
+summary(pb)
+plot(pb)
+
+system.time(
+ pm <- probe.match(
+ pb,
+ est=c("log.r","log.phi","N.0"),
+ parscale=c(0.1,0.1,0.1),
+ nsim=1000,
+ seed=838775L,
+ method="Nelder-Mead",
+ reltol=1e-7,
+ fail.value=1e9
+ )
+ )
+plot(pm)
+
+cbind(truth=coef(ricker),est=coef(pm),guess=coef(po))
+
+pb <- probe(
+ po,
probes=probe.marginal(
var="y",
transform=sqrt,
@@ -143,9 +202,15 @@
probe.acf(
var="y",
transform=sqrt,
- lag.max=0,
+ lag.max=2,
type="cov"
- )
+ ),
+ probe.nlar(
+ var="y",
+ transform=sqrt,
+ lags=c(1,2),
+ powers=1
+ )
),
nsim=1000,
seed=838775L
@@ -194,5 +259,4 @@
)
)
-
dev.off()
More information about the pomp-commits
mailing list