[Pomp-commits] r324 - in pkg: . R inst inst/doc man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 24 17:58:43 CEST 2010
Author: kingaa
Date: 2010-09-24 17:58:43 +0200 (Fri, 24 Sep 2010)
New Revision: 324
Added:
pkg/src/Makevars
pkg/src/mat.c
pkg/src/probe_acf.c
pkg/src/probe_marginal.c
pkg/tests/ricker-probe.R
pkg/tests/ricker-probe.Rout.save
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/pomp.bib
pkg/man/basic-probes.Rd
pkg/src/pomp_internal.h
pkg/tests/ou2-probe-match.R
pkg/tests/ou2-probe-match.Rout.save
Log:
- probe.marginal and probe.median are introduced, based on the recommendation of Wood (2010)
- probe.acf and probe.marginal are now coded in C for speed
- probe.acf is slightly changed to use the natural definition of ACF
- these ideas follow the recommendations of Simon Wood (2010)
- a test file, using the Ricker example, is included
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-09-24 13:57:08 UTC (rev 323)
+++ pkg/DESCRIPTION 2010-09-24 15:58:43 UTC (rev 324)
@@ -1,7 +1,7 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.32-7
+Version: 0.33-1
Date: 2010-09-24
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing,
Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-09-24 13:57:08 UTC (rev 323)
+++ pkg/NAMESPACE 2010-09-24 15:58:43 UTC (rev 324)
@@ -12,6 +12,9 @@
pfilter_computations,
apply_probe_data,
apply_probe_sim,
+ probe_marginal_setup,
+ probe_marginal_solve,
+ probe_acf,
do_rprocess,
do_dprocess,
do_rmeasure,
@@ -61,6 +64,7 @@
nlf,
traj.match,
probe.mean,
+ probe.median,
probe.var,
probe.sd,
probe.period,
@@ -68,6 +72,7 @@
probe.acf,
probe.cov,
probe.cor,
+ probe.marginal,
probe.match,
spect.match
)
Modified: pkg/R/basic-probes.R
===================================================================
--- pkg/R/basic-probes.R 2010-09-24 13:57:08 UTC (rev 323)
+++ pkg/R/basic-probes.R 2010-09-24 15:58:43 UTC (rev 324)
@@ -1,16 +1,29 @@
probe.mean <- function (var, trim = 0, transform = identity, na.rm = TRUE) {
+ if (length(var)>1) stop(sQuote("probe.mean")," is a univariate probe")
+ transform <- match.fun(transform)
function(y) mean(x=transform(y[var,]),trim=trim,na.rm=na.rm)
}
+probe.median <- function (var, na.rm = TRUE) {
+ if (length(var)>1) stop(sQuote("probe.median")," is a univariate probe")
+ function(y) median(x=as.numeric(y[var,]),na.rm=na.rm)
+}
+
probe.var <- function (var, transform = identity, na.rm = TRUE) {
+ if (length(var)>1) stop(sQuote("probe.var")," is a univariate probe")
+ transform <- match.fun(transform)
function(y) var(x=transform(y[var,]),na.rm=na.rm)
}
probe.sd <- function (var, transform = identity, na.rm = TRUE) {
+ if (length(var)>1) stop(sQuote("probe.sd")," is a univariate probe")
+ transform <- match.fun(transform)
function(y) sd(x=transform(y[var,]),na.rm=na.rm)
}
probe.period <- function (var, kernel.width, transform = identity) {
+ if (length(var)>1) stop(sQuote("probe.period")," is a univariate probe")
+ transform <- match.fun(transform)
function (y) {
zz <- spec.pgram(
x=transform(y[var,]),
@@ -26,29 +39,11 @@
}
probe.quantile <- function (var, prob, transform = identity) {
+ if (length(var)>1) stop(sQuote("probe.quantile")," is a univariate probe")
+ transform <- match.fun(transform)
function (y) quantile(transform(y[var,]),probs=prob)
}
-probe.acf <- function (
- var,
- lag = 0,
- type = c("correlation", "covariance", "partial"),
- transform = identity,
- ...
- ) {
- args <- list(...)
- type <- match.arg(type)
- mlag <- max(lag)
- function (y) {
- zz <- do.call(acf,c(list(x=transform(y[var,]),lag.max=mlag,plot=FALSE,type=type),args))
- if (type=="partial")
- val <- zz$acf[lag]
- else
- val <- zz$acf[lag+1]
- val
- }
-}
-
probe.cov <- function (
vars,
lag,
@@ -57,6 +52,7 @@
) {
method <- match.arg(method)
lag <- as.integer(lag)
+ transform <- match.fun(transform)
var1 <- vars[1]
if (length(vars)>1)
var2 <- vars[2]
@@ -88,6 +84,7 @@
) {
method <- match.arg(method)
lag <- as.integer(lag)
+ transform <- match.fun(transform)
var1 <- vars[1]
if (length(vars)>1)
var2 <- vars[2]
@@ -110,3 +107,27 @@
val
}
}
+
+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)
+ setup <- .Call(probe_marginal_setup,transform(ref),order,diff)
+ function (y) .Call(
+ probe_marginal_solve,
+ x=transform(y[var,]),
+ setup=setup,
+ diff=diff
+ )
+}
+
+probe.acf <- function (var, lag.max, type = c("covariance", "correlation"), transform = identity) {
+ type <- match.arg(type)
+ transform <- match.fun(transform)
+ corr <- type=="correlation"
+ function (y) .Call(
+ probe_acf,
+ x=transform(y[var,,drop=FALSE]),
+ lag_max=lag.max,
+ corr=corr
+ )
+}
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-09-24 13:57:08 UTC (rev 323)
+++ pkg/inst/ChangeLog 2010-09-24 15:58:43 UTC (rev 324)
@@ -1,5 +1,23 @@
+2010-09-24 kingaa
+
+ * [r321] R/pomp-methods.R, man/pomp-methods.Rd: - add S3 method
+ as.data.frame.pomp
+ * [r320] DESCRIPTION, src/pfilter.c: - bug fix
+
+2010-09-23 kingaa
+
+ * [r319] R/probe.R: - empirical p-values are computed incorrectly.
+ fix this.
+
+2010-09-21 kingaa
+
+ * [r317] inst/doc/advanced_topics_in_pomp.pdf,
+ inst/doc/intro_to_pomp.pdf, inst/doc/ou2-first-mif.rda,
+ inst/doc/ou2-trajmatch.rda: - update the vignettes
+
2010-09-20 kingaa
+ * [r316] inst/ChangeLog: - update changelog
* [r315] DESCRIPTION, NAMESPACE, R/pfilter.R, inst/ChangeLog,
inst/doc/advanced_topics_in_pomp.pdf, inst/doc/intro_to_pomp.pdf,
src/pfilter.c, tests/ricker.Rout.save: - move most time-consuming
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/pomp.bib
===================================================================
--- pkg/inst/doc/pomp.bib 2010-09-24 13:57:08 UTC (rev 323)
+++ pkg/inst/doc/pomp.bib 2010-09-24 15:58:43 UTC (rev 324)
@@ -3,7 +3,7 @@
@ARTICLE{Andrieu2010,
author = {Andrieu, Christophe and Doucet, Arnaud and Holenstein, Roman},
- title = {Particle {M}arkov chain {M}onte {C}arlo methods},
+ title = {{P}article {M}arkov chain {M}onte {C}arlo methods},
journal = {Journal of the Royal Statistical Society, Series B},
year = {2010},
volume = {72},
@@ -31,8 +31,8 @@
@ARTICLE{Arulampalam2002,
author = {Arulampalam, M. S. and Maskell, S. and Gordon, N. and Clapp, T.},
- title = {A Tutorial on Particle Filters for Online Nonlinear, Non-{G}aussian
- {B}ayesian Tracking},
+ title = {{A} {T}utorial on {P}article {F}ilters for {O}nline {N}onlinear,
+ {N}on-{G}aussian {B}ayesian {T}racking},
journal = {IEEE Transactions on Signal Processing},
year = {2002},
volume = {50},
@@ -45,8 +45,8 @@
@ARTICLE{Bhadra2010,
author = {Bhadra, Anindya},
- title = {Discussion of `Particle {M}arkov chain {M}onte {C}arlo methods' by
- C.\ Andrieu, A.\ Doucet and R.\ Holenstein},
+ title = {{D}iscussion of `{P}article {M}arkov chain {M}onte {C}arlo methods'
+ by {C}.\ {A}ndrieu, {A}.\ {D}oucet and {R}.\ {H}olenstein},
journal = {Journal of the Royal Statistical Society, Series B},
year = {2010},
volume = {72},
@@ -60,7 +60,7 @@
@ARTICLE{Breto2009,
author = {Carles Bret\'{o} and Daihai He and Edward L. Ionides and Aaron A.
King},
- title = {Time series analysis via mechanistic models},
+ title = {{T}ime series analysis via mechanistic models},
journal = {Annals of Applied Statistics},
year = {2009},
volume = {3},
@@ -92,7 +92,7 @@
@ARTICLE{Ellner1998,
author = {S. P. Ellner and B. A. Bailey and G. V. Bobashev and A. R. Gallant
and B. T. Grenfell and D. W. Nychka},
- title = {Noise and nonlinearity in measles epidemics: Combining mechanistic
+ title = {{N}oise and nonlinearity in measles epidemics: {C}ombining mechanistic
and statistical approaches to population modeling},
journal = {American Naturalist},
year = {1998},
@@ -128,7 +128,7 @@
@ARTICLE{Gillespie1977a,
author = {D. T. Gillespie},
- title = {Exact Stochastic Simulation of Coupled Chemical Reactions},
+ title = {{E}xact {S}tochastic {S}imulation of {C}oupled {C}hemical {R}eactions},
journal = {Journal of Physical Chemistry},
year = {1977},
volume = {81},
@@ -140,8 +140,8 @@
@ARTICLE{He2010,
author = {He, Daihai and Ionides, Edward L. and King, Aaron A.},
- title = {Plug-and-play inference for disease dynamics: measles in large and
- small populations as a case study},
+ title = {{P}lug-and-play inference for disease dynamics: measles in large
+ and small populations as a case study},
journal = {Journal of the Royal Society Interface},
year = {2010},
volume = {7},
@@ -173,7 +173,7 @@
@ARTICLE{Ionides2006,
author = {E. L. Ionides and C. Bret{\'o} and Aaron A. King},
- title = {Inference for nonlinear dynamical systems},
+ title = {{I}nference for nonlinear dynamical systems},
journal = {Proceedings of the National Academy of Sciences of the U.S.A.},
year = {2006},
volume = {103},
@@ -204,7 +204,7 @@
@ARTICLE{Kendall1999,
author = {B. E. Kendall and C. J. Briggs and W. W. Murdoch and P. Turchin and
S. P. Ellner and E. McCauley and R. M. Nisbet and S. N. Wood},
- title = {Why do populations cycle? A synthesis of statistical and mechanistic
+ title = {{W}hy do populations cycle? {A} synthesis of statistical and mechanistic
modeling approaches},
journal = {Ecology},
year = {1999},
@@ -248,8 +248,8 @@
@ARTICLE{Kendall2005,
author = {Kendall, B. E. and Ellner, S. P. and McCauley, E. and Wood, S. N.
and Briggs, C. J. and Murdoch, W. M. and Turchin, P.},
- title = {Population cycles in the pine looper moth: Dynamical tests of mechanistic
- hypotheses},
+ title = {{P}opulation cycles in the pine looper moth: {D}ynamical tests of
+ mechanistic hypotheses},
journal = {Ecological Monographs},
year = {2005},
volume = {75},
@@ -297,7 +297,7 @@
@ARTICLE{King2008,
author = {King, Aaron A. and Ionides, Edward L. and Pascual, Mercedes and Bouma,
Menno J.},
- title = {Inapparent infections and cholera dynamics},
+ title = {{I}napparent infections and cholera dynamics},
journal = {Nature},
year = {2008},
volume = {454},
@@ -340,7 +340,8 @@
@INCOLLECTION{Liu2001b,
author = {Liu, J and West, M.},
- title = {Combining Parameter and State Estimation in Simulation-Based Filtering},
+ title = {{C}ombining {P}arameter and {S}tate {E}stimation in {S}imulation-{B}ased
+ {F}iltering},
booktitle = {Sequential {M}onte {C}arlo Methods in Practice},
publisher = {Springer, New York},
year = {2001},
@@ -353,7 +354,7 @@
@ARTICLE{Reuman2006,
author = {Reuman, D. C. and Desharnais, R. A. and Costantino, R. F. and Ahmad,
O. S. and Cohen, J. E.},
- title = {Power spectra reveal the influence of stochasticity on nonlinear
+ title = {{P}ower spectra reveal the influence of stochasticity on nonlinear
population dynamics},
journal = {Proceedings of the National Academy of Sciences of the U.S.A.},
year = {2006},
@@ -364,7 +365,7 @@
@ARTICLE{Wearing2009,
author = {Helen J Wearing and Pejman Rohani},
- title = {Estimating the duration of pertussis immunity using epidemiological
+ title = {{E}stimating the duration of pertussis immunity using epidemiological
signatures.},
journal = {PLoS Pathogens},
year = {2009},
@@ -398,6 +399,21 @@
timestamp = {2010.02.05}
}
+ at ARTICLE{Wood2010,
+ author = {Wood, Simon N.},
+ title = {{S}tatistical inference for noisy nonlinear ecological dynamic systems},
+ journal = {Nature},
+ year = {2010},
+ volume = {466},
+ pages = {1102--1104},
+ month = aug,
+ doi = {10.1038/nature09319},
+ file = {Wood2010.pdf:Wood2010.pdf:PDF;Wood2010_Supplement.pdf:Wood2010_Supplement.pdf:PDF},
+ issn = {1476-4687},
+ owner = {kingaa},
+ timestamp = {2010.08.17}
+}
+
@comment{jabref-meta: selector_publisher:}
@comment{jabref-meta: selector_author:}
Modified: pkg/man/basic-probes.Rd
===================================================================
--- pkg/man/basic-probes.Rd 2010-09-24 13:57:08 UTC (rev 323)
+++ pkg/man/basic-probes.Rd 2010-09-24 15:58:43 UTC (rev 324)
@@ -1,6 +1,7 @@
\name{basic.probes}
\alias{basic.probes}
\alias{probe.mean}
+\alias{probe.median}
\alias{probe.var}
\alias{probe.sd}
\alias{probe.period}
@@ -8,6 +9,7 @@
\alias{probe.acf}
\alias{probe.cov}
\alias{probe.cor}
+\alias{probe.marginal}
\title{Some probes for partially-observed Markov processes}
\description{
Several simple and configurable probes are provided in the package.
@@ -15,20 +17,22 @@
}
\usage{
probe.mean(var, trim = 0, transform = identity, na.rm = TRUE)
+probe.median(var, na.rm = TRUE)
probe.var(var, transform = identity, na.rm = TRUE)
probe.sd(var, transform = identity, na.rm = TRUE)
-probe.period(var, kernel.width, transform = identity)
-probe.quantile(var, prob, transform = identity)
-probe.acf(var, lag = 0, type = c("correlation", "covariance", "partial"),
- transform = identity, \dots)
+probe.marginal(var, ref, order = 3, diff = 1, transform = identity)
+probe.acf(var, lag.max, type = c("covariance", "correlation"),
+ transform = identity)
probe.cov(vars, lag, method = c("pearson", "kendall", "spearman"),
transform = identity)
probe.cor(vars, lag, method = c("pearson", "kendall", "spearman"),
transform = identity)
+probe.period(var, kernel.width, transform = identity)
+probe.quantile(var, prob, transform = identity)
}
\arguments{
\item{var}{
- character; the name or index of the observed variable.
+ character; the name(s) of the observed variable(s).
}
\item{trim}{
the fraction of observations to be trimmed (see \code{\link{mean}}).
@@ -46,19 +50,34 @@
\item{prob}{
a single probability; the quantile to compute: see \code{\link{quantile}}.
}
+ \item{lag.max}{
+ The maximum lag at which the ACF is computed.
+ }
\item{lag}{
- In \code{probe.acf}, the lag at which the ACF is computed.
In \code{probe.cov} and \code{probe.cor}, the lag between time series.
}
\item{type}{
- The type of ACF to compute: see \code{\link{acf}}.
+ Compute autocorrelation or autocovariance?
}
\item{vars}{
- length-2 character vector; the names or indices of the observed variables between which covariance or correlation is to be computed.
+ length-2 character vector;
+ the names or indices of the observed variables between which covariance or correlation is to be computed.
}
\item{method}{
The type of covariance or correlation to compute: see \code{\link{cov}} and \code{\link{cor}}.
}
+ \item{ref}{
+ empirical reference distribution.
+ Simulated data will be regressed against the values of \code{ref}, sorted and, optionally, differenced.
+ The resulting regression coefficients capture information about the shape of the marginal distribution.
+ A good choice for \code{ref} is the data itself.
+ }
+ \item{order}{
+ order of polynomial regression.
+ }
+ \item{diff}{
+ order of differencing to perform.
+ }
\item{\dots}{
Additional arguments to be passed through to the probe computation.
}
@@ -71,26 +90,39 @@
Each of these functions is relatively simple.
See the source code for a complete understanding of what each does.
\describe{
- \item{\code{probe.mean}, \code{probe.var}, \code{probe.sd}}{
- compute the mean, variance, and standard deviation of variable \code{var}, respectively.
+ \item{\code{probe.mean}, \code{probe.median}, \code{probe.var}, \code{probe.sd}}{
+ return functions that compute the mean, median, variance, and standard deviation of variable \code{var}, respectively.
}
\item{\code{probe.period}}{
- estimates the period of the Fourier component of the \code{var} series with largest power.
+ returns a function that estimates the period of the Fourier component of the \code{var} series with largest power.
}
- \item{\code{probe.quantile}}{
- estimates the \code{prob}-th quantile of variable \code{var}.
+ \item{\code{probe.marginal}}{
+ returns a function that
+ regresses the marginal distribution of variable \code{var} against the reference distribution \code{ref}.
+ If \code{diff>0}, the data and the reference distribution are first differenced \code{diff} times and centered.
+ Polynomial regression of order \code{order} is used.
+ This probe returns \code{order} regression coefficients (the intercept is zero).
}
\item{\code{probe.acf}}{
- computes the autocorrelation or autocovariance function of variable \code{var} at lag \code{lag}.
+ returns a function that,
+ if \code{type=="covariance"}, computes the autocovariance function of variable \code{var} at lags 0 through \code{lag.max};
+ if \code{type=="correlation"}, computes the autocorrelation function of variable \code{var} at lags 1 through \code{lag.max}.
}
\item{\code{probe.cov}, \code{probe.cor}}{
- compute the covariance and correlation, respectively, between the two variables named in \code{vars} at relative lag \code{lag}.}
+ return functions that compute the covariance and correlation, respectively, between the two variables named in \code{vars} at relative lag \code{lag}.}
+ \item{\code{probe.quantile}}{
+ returns a function that estimates the \code{prob}-th quantile of variable \code{var}.
+ }
}
}
\references{
B. E. Kendall, C. J. Briggs, W. M. Murdoch, P. Turchin, S. P. Ellner, E. McCauley, R. M. Nisbet, S. N. Wood
Why do populations cycle? A synthesis of statistical and mechanistic modeling approaches,
Ecology, 80:1789--1805, 1999.
+
+ S. N. Wood
+ Statistical inference for noisy nonlinear ecological dynamic systems,
+ Nature, 466: 1102--1104, 2010.
}
\author{
Daniel C. Reuman (d.reuman at imperial dot ac dot uk)
Added: pkg/src/Makevars
===================================================================
--- pkg/src/Makevars (rev 0)
+++ pkg/src/Makevars 2010-09-24 15:58:43 UTC (rev 324)
@@ -0,0 +1 @@
+ PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
Property changes on: pkg/src/Makevars
___________________________________________________________________
Added: svn:executable
+ *
Added: pkg/src/mat.c
===================================================================
--- pkg/src/mat.c (rev 0)
+++ pkg/src/mat.c 2010-09-24 15:58:43 UTC (rev 324)
@@ -0,0 +1,60 @@
+// -*- mode: C++; -*-
+// some functions for linear algebra written by Simon N. Wood.
+// codes have been modified by AAK to suit his own peccadilloes
+
+#include <math.h>
+#include <R.h>
+#include <R_ext/Linpack.h>
+#include <R_ext/Lapack.h>
+#include "pomp_internal.h"
+
+void pomp_backsolve (double *R, int *r, int *c, double *B, double *C, int *bc) {
+ double x;
+ int i, j, k;
+ for (j = 0; j < *bc; j++) { // work across columns of B & C
+ for (i = *c-1; i >= 0; i--) { // work up each column of B & C
+ for (k = i+1, x = 0.0; k < *c; k++)
+ x += R[i+*r*k]*C[k+*c*j];
+ C[i+*c*j] = (B[i+*c*j]-x)/R[i+*r*i];
+ }
+ }
+}
+
+void pomp_qr (double *x, int *r, int *c, int *pivot, double *tau) {
+ int info, lwork = -1, i;
+ double work1, *work;
+ // workspace query
+ F77_NAME(dgeqp3)(r,c,x,r,pivot,tau,&work1,&lwork,&info);
+ lwork = (int) floor(work1);
+ if ((work1-lwork) >0.5) lwork++;
+ work = (double *) Calloc(lwork,double);
+ // actual call
+ F77_NAME(dgeqp3)(r,c,x,r,pivot,tau,work,&lwork,&info);
+ Free(work);
+ for (i = 0; i < *c; i++) pivot[i]--;
+ // ... for 'tis C in which we work and not the 'cursed Fortran...
+}
+
+
+void pomp_qrqy (double *b, double *a, double *tau, int *r, int *c, int *k, int *left, int *tp) {
+ char side, trans='N';
+ int lda, lwork = -1, info;
+ double *work, work1;
+
+ if (! *left) {
+ side='R';
+ lda = *c;
+ } else {
+ side = 'L';
+ lda= *r;
+ }
+ if (*tp) trans='T';
+ // workspace query
+ F77_NAME(dormqr)(&side,&trans,r,c,k,a,&lda,tau,b,r,&work1,&lwork,&info);
+ lwork = (int) floor(work1);
+ if ((work1-lwork) > 0.5) lwork++;
+ work = (double *) Calloc(lwork,double);
+ // actual call
+ F77_NAME(dormqr)(&side,&trans,r,c,k,a,&lda,tau,b,r,work,&lwork,&info);
+ Free(work);
+}
Property changes on: pkg/src/mat.c
___________________________________________________________________
Added: svn:eol-style
+ native
Modified: pkg/src/pomp_internal.h
===================================================================
--- pkg/src/pomp_internal.h 2010-09-24 13:57:08 UTC (rev 323)
+++ pkg/src/pomp_internal.h 2010-09-24 15:58:43 UTC (rev 324)
@@ -149,4 +149,8 @@
#endif
+void pomp_backsolve(double *, int *, int *, double *, double *, int *);
+void pomp_qr(double *, int *, int *, int *, double *);
+void pomp_qrqy(double *, double *, double *, int *, int *, int *, int *, int *);
+
#endif
Added: pkg/src/probe_acf.c
===================================================================
--- pkg/src/probe_acf.c (rev 0)
+++ pkg/src/probe_acf.c 2010-09-24 15:58:43 UTC (rev 324)
@@ -0,0 +1,98 @@
+// -*- mode: C++ -*-
+
+#include "pomp_internal.h"
+#include <stdio.h>
+
+static void pomp_acf (double *acf, double *x, int n, int nvars, int maxlag, int correlation);
+
+SEXP probe_acf (SEXP x, SEXP lag_max, SEXP corr) {
+ int nprotect = 0;
+ SEXP acf, acf_names, cacf;
+ SEXP X_names, Dim;
+ int maxlag, correlation, nvars, n;
+ int j, k, l;
+ double *p, *p1;
+ char tmp[BUFSIZ];
+
+ maxlag = *(INTEGER(AS_INTEGER(lag_max))); // maximum lag
+ correlation = *(INTEGER(AS_INTEGER(corr))); // correlation, or covariance?
+
+ PROTECT(Dim = GET_DIM(x)); nprotect++;
+ nvars = INTEGER(Dim)[0]; // nvars = # of variables
+ n = INTEGER(Dim)[1]; // n = # of observations
+ // depending on how this function is called, it may be necessary to duplicate x in the next line
+ PROTECT(x = AS_NUMERIC(x)); nprotect++;
+ PROTECT(X_names = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
+
+ PROTECT(acf = NEW_NUMERIC((maxlag+1)*nvars)); nprotect++;
+
+ pomp_acf(REAL(acf),REAL(x),n,nvars,maxlag,correlation);
+
+ if (correlation) {
+ PROTECT(cacf = NEW_NUMERIC(maxlag*nvars)); nprotect++;
+ PROTECT(acf_names = NEW_STRING(maxlag*nvars)); nprotect++;
+ for (j = 0, l = 0, p = REAL(cacf), p1 = REAL(acf)+1; j < nvars; j++, p1++) {
+ for (k = 1; k <= maxlag; k++, p++, p1++, l++) {
+ *p = *p1; // copy correlations from acf into cacf
+ snprintf(tmp,BUFSIZ,"acf.%ld.%s",k,CHARACTER_DATA(STRING_ELT(X_names,j)));
+ SET_STRING_ELT(acf_names,l,mkChar(tmp));
+ }
+ }
+ SET_NAMES(cacf,acf_names);
+ } else {
+ PROTECT(acf_names = NEW_STRING((maxlag+1)*nvars)); nprotect++;
+ for (j = 0, l = 0; j < nvars; j++) {
+ for (k = 0; k <= maxlag; k++, l++) {
+ snprintf(tmp,BUFSIZ,"acf.%ld.%s",k,CHARACTER_DATA(STRING_ELT(X_names,j)));
+ SET_STRING_ELT(acf_names,l,mkChar(tmp));
+ }
+ }
+ SET_NAMES(acf,acf_names);
+ }
+
+ UNPROTECT(nprotect);
+ if (correlation) return(cacf);
+ else return(acf);
+}
+
+// vectorized routine for ACF calculation
+// thanks to Simon N. Wood for the original version of this code
+// modifications due to AAK
+static void pomp_acf (double *acf, double *x, int n, int nvars, int maxlag, int correlation) {
+ double sum, *p, *p0, *p1, *p2;
+ int j, k, lag, ct;
+
+ // first center each row
+ for (j = 0, p = x; j < nvars; j++, p++) {
+ for (k = 0, p0 = p, sum = 0, ct = 0; k < n; p0 += nvars, k++) {
+ if (R_FINITE(*p0)) {
+ sum += *p0;
+ ct++;
+ }
+ }
+ if (ct < 1) error("series %ld has no data",j);
+ sum /= ((double) ct); // mean of x[j,]
+ for (k = 0, p0 = p; k < n; p0 += nvars, k++)
+ if (R_FINITE(*p0)) *p0 -= sum;
+ }
+
+ // compute covariances
+ for (j = 0, p0 = x, p = acf; j < nvars; j++, p0++) { // loop over series
+ for (lag = 0; lag <= maxlag; lag++, p++) { // loop over lags
+ for (k = 0, ct = 0, sum = 0, p1 = p0, p2 = p0+lag*nvars; k < n-lag; k++, p1 += nvars, p2 += nvars)
+ if (R_FINITE(*p1) && R_FINITE(*p2)) {
+ sum += (*p1)*(*p2);
+ ct++;
+ }
+ *p = (ct > 0) ? sum/((double) ct) : R_NaReal;
+ //// *p = (ct > 0) ? sum/((double) n) : R_NaReal; // strangely, this is apparently what R's 'acf' function does
+ }
+ }
+
+ // convert to correlations if desired
+ if (correlation)
+ for (j = 0, p = acf; j < nvars; j++, p += maxlag+1)
+ for (lag = 0, p0 = p, sum = *p; lag <= maxlag; lag++, p0++)
+ *p0 /= sum;
+
+}
Property changes on: pkg/src/probe_acf.c
___________________________________________________________________
Added: svn:eol-style
+ native
Added: pkg/src/probe_marginal.c
===================================================================
--- pkg/src/probe_marginal.c (rev 0)
+++ pkg/src/probe_marginal.c 2010-09-24 15:58:43 UTC (rev 324)
@@ -0,0 +1,148 @@
+// -*- mode: C++ -*-
+
+#include "pomp_internal.h"
+#include <stdio.h>
+#include <R_ext/Linpack.h>
+#include <R_ext/Lapack.h>
+
+static void order_reg_solve (double *beta, double *x, double *mm, double *tau, int *pivot, int n, int np, int diff);
+static void order_reg_model_matrix (double *z, double *X, double *pivot, double *tau, int n, int np, int diff);
+
+SEXP probe_marginal_setup (SEXP ref, SEXP order, SEXP diff) {
+ int nprotect = 0;
+ SEXP z, mm, tau, pivot, retval, retvalnames;
+ int n, nx, np, df, dim[2];
+
+ np = *(INTEGER(AS_INTEGER(order))); // order of polynomial regression
+ df = *(INTEGER(AS_INTEGER(diff))); // order of differencing
+ n = LENGTH(ref); // n = number of observations
+ nx = n - df; // nx = rows in model matrix
+ dim[0] = nx; dim[1] = np; // dimensions of model matrix
+
+ if (nx < 1) error("must have diff < number of observations");
+
+ PROTECT(z = duplicate(AS_NUMERIC(ref))); nprotect++;
+ PROTECT(mm = makearray(2,dim)); nprotect++;
+ PROTECT(tau = NEW_NUMERIC(np)); nprotect++;
+ PROTECT(pivot = NEW_INTEGER(np)); nprotect++;
+
+ PROTECT(retval = NEW_LIST(3)); nprotect++;
+ PROTECT(retvalnames = NEW_CHARACTER(3)); nprotect++;
+ SET_STRING_ELT(retvalnames,0,mkChar("mm"));
+ SET_STRING_ELT(retvalnames,1,mkChar("tau"));
+ SET_STRING_ELT(retvalnames,2,mkChar("pivot"));
+ SET_ELEMENT(retval,0,mm);
+ SET_ELEMENT(retval,1,tau);
+ SET_ELEMENT(retval,2,pivot);
+ SET_NAMES(retval,retvalnames);
+
+ order_reg_model_matrix(REAL(z),REAL(mm),REAL(tau),INTEGER(pivot),n,np,df);
+
+ UNPROTECT(nprotect);
+ return(retval);
+}
+
+SEXP probe_marginal_solve (SEXP x, SEXP setup, SEXP diff) {
+ int nprotect = 0;
+ SEXP dimx, dimm;
+ SEXP X, mm, tau, pivot, beta, beta_names;
+ int n, nx, np, df;
+ int i;
+ char tmp[BUFSIZ];
+
+ df = *(INTEGER(AS_INTEGER(diff))); // order of differencing
+ n = LENGTH(x); // n = number of observations
+
+ // unpack the setup information
+ PROTECT(mm = VECTOR_ELT(setup,0)); nprotect++; // QR-decomposed model matrix
+ PROTECT(tau = VECTOR_ELT(setup,1)); nprotect++; // diagonals
+ PROTECT(pivot = VECTOR_ELT(setup,2)); nprotect++; // pivots
+
+ PROTECT(dimm = GET_DIM(mm)); nprotect++;
+ nx = INTEGER(dimm)[0]; // nx = number of rows in model matrix
+ np = INTEGER(dimm)[1]; // np = order of polynomial
+
+ if (n-df != nx) error("length of 'ref' must equal length of data");
+ PROTECT(X = duplicate(AS_NUMERIC(x))); nprotect++;
+
+ PROTECT(beta = NEW_NUMERIC(np)); nprotect++;
+ PROTECT(beta_names = NEW_STRING(np)); nprotect++;
+ for (i = 0; i < np; i++) {
+ snprintf(tmp,BUFSIZ,"od.%ld",i+1);
+ SET_STRING_ELT(beta_names,i,mkChar(tmp));
+ }
+ SET_NAMES(beta,beta_names);
+
+ order_reg_solve(REAL(beta),REAL(X),REAL(mm),REAL(tau),INTEGER(pivot),n,np,df);
+
+ UNPROTECT(nprotect);
+ return(beta);
+}
+
+// thanks to Simon N. Wood for the original version of the following code
+static void order_reg_model_matrix (double *z, double *X, double *tau, double *pivot, int n, int np, int diff) {
+ // z is an n vector, containing no NAs
+ // X is an (n-diff) by np double matrix
+ // pivot is an (n-diff) integer vector
+ // tau is an (n-diff) double vector
+ // This routine first differences z 'diff' times.
+ // z is centred and sorted
+ // then the model matrix is set up and QR decomposed
+ int nx;
+ double *X1, *X2, xx;
+ int i, j;
+
+ for (i = 0, nx = n; i < diff; i++) { // differencing loop
+ nx--;
+ for (j = 0; j < nx; j++) z[j] = z[j+1] - z[j];
+ }
+ // nx = number of rows in model matrix
+ // z is now an nx-vector
+
+ // center z
+ for (j = 0, xx = 0.0; j < nx; j++) xx += z[j]; xx /= nx; // xx = mean(z)
+ for (j = 0; j < nx; j++) z[j] -= xx;
+
+ // now sort
+ R_qsort(z,1,nx);
+
+ // Now create the model matrix, using contents of z
+ if (np < 1) np = 1;
+ for (j = 0; j < nx; j++) X[j] = z[j]; // first column
+ for (i = 1, X1 = X, X2 = X+nx; i < np; i++, X1 += nx, X2 += nx)
+ for (j = 0; j < nx; j++) X2[j] = X1[j]*z[j];
+
+ // QR decompose the model matrix
+ for (i = 0; i < np; i++) pivot[i] = 0;
+ pomp_qr(X,&nx,&np,pivot,tau);
+
+}
+
+// thanks to Simon N. Wood for the original version of the following code
+static void order_reg_solve (double *beta, double *x, double *mm, double *tau, int *pivot, int n, int np, int diff) {
+ int nx, one = 1;
+ double xx, coef[np];
+ int i, j;
+
+ for (i = 0, nx = n; i < diff; i++) { // differencing loop
+ nx--;
+ for (j = 0; j < nx; j++) x[j] = x[j+1] - x[j];
+ }
+ // nx = number of rows in model matrix
+ // x is now an nx-vector
+
+ // center x
+ for (j = 0, xx = 0.0; j < nx; j++) xx += x[j]; xx /= nx; // xx = mean(x)
+ for (j = 0; j < nx; j++) x[j] -= xx;
+
+ // now sort
+ R_qsort(x,1,nx);
+
+ // solve R b = Q'x for b
+ pomp_qrqy(x,mm,tau,&nx,&one,&np,&one,&one); // y <- Q'y
+ pomp_backsolve(mm,&nx,&np,x,coef,&one); // b <- R^{-1} Q'y
+
+ // unpivot beta
+ for (i = 0; i < np; i++) beta[i] = coef[pivot[i]];
+
+}
Property changes on: pkg/src/probe_marginal.c
___________________________________________________________________
Added: svn:eol-style
+ native
Modified: pkg/tests/ou2-probe-match.R
===================================================================
--- pkg/tests/ou2-probe-match.R 2010-09-24 13:57:08 UTC (rev 323)
+++ pkg/tests/ou2-probe-match.R 2010-09-24 15:58:43 UTC (rev 324)
@@ -35,7 +35,7 @@
pm.ou2 <- probe(
ou2,
probes=list(
- y1.lag3=probe.acf(var="y1",lag=3,type="corr"),
+ y1acf=probe.acf(var="y1",lag.max=2,type="corr"),
y2.cov12=probe.cov(vars=c("y1"),lag=12,method="spearman"),
y12.cov8=probe.cov(vars=c("y2","y1"),lag=8,method="pearson")
),
@@ -48,7 +48,7 @@
ou2,
probes=list(
y1=probe.quantile(var="y1",prob=seq(0.1,0.9,by=0.1)),
- acf=probe.acf(var="y2",lag=c(2,4,7),transform=identity),
+ probe.acf(var=c("y1","y2"),lag.max=4,transform=identity),
pd=probe.period(var="y1",kernel.width=3)
),
nsim=200
Modified: pkg/tests/ou2-probe-match.Rout.save
===================================================================
--- pkg/tests/ou2-probe-match.Rout.save 2010-09-24 13:57:08 UTC (rev 323)
+++ pkg/tests/ou2-probe-match.Rout.save 2010-09-24 15:58:43 UTC (rev 324)
@@ -1,8 +1,7 @@
-R version 2.12.0 Under development (unstable) (2010-08-18 r52775)
+R version 2.11.1 (2010-05-31)
Copyright (C) 2010 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.
@@ -58,7 +57,7 @@
$pvals
y1.mean y2.mean y1.sd y2.sd
-0.20359281 0.93812375 0.05189621 0.10778443
+0.20359281 0.94211577 0.05189621 0.10778443
> summary(pm.po)
$coef
@@ -73,8 +72,8 @@
0.054 0.644 1.000 1.000
$pvals
- y1.mean y2.mean y1.sd y2.sd
-0.1117764 0.7105788 0.0000000 0.0000000
+ y1.mean y2.mean y1.sd y2.sd
+0.111776447 0.714570858 0.003992016 0.003992016
>
> plot(pm.ou2)
@@ -83,7 +82,7 @@
> pm.ou2 <- probe(
+ ou2,
+ probes=list(
-+ y1.lag3=probe.acf(var="y1",lag=3,type="corr"),
++ y1acf=probe.acf(var="y1",lag.max=2,type="corr"),
+ y2.cov12=probe.cov(vars=c("y1"),lag=12,method="spearman"),
+ y12.cov8=probe.cov(vars=c("y2","y1"),lag=8,method="pearson")
+ ),
@@ -99,19 +98,19 @@
[1] 500
$quantiles
- y1.lag3 y2.cov12 y12.cov8
- 0.110 0.726 0.462
+y1acf.acf.1.y1 y1acf.acf.2.y1 y2.cov12 y12.cov8
+ 0.092 0.054 0.726 0.462
$pvals
- y1.lag3 y2.cov12 y12.cov8
-0.2235529 0.5469062 0.9261477
+y1acf.acf.1.y1 y1acf.acf.2.y1 y2.cov12 y12.cov8
+ 0.1876248 0.1117764 0.5508982 0.9261477
>
> pb <- probe(
+ ou2,
+ probes=list(
+ y1=probe.quantile(var="y1",prob=seq(0.1,0.9,by=0.1)),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 324
More information about the pomp-commits
mailing list