[Pomp-commits] r387 - in pkg: . R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 19 15:18:09 CEST 2010
Author: kingaa
Date: 2010-10-19 15:18:08 +0200 (Tue, 19 Oct 2010)
New Revision: 387
Added:
pkg/src/synth_lik.c
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/basic-probes.R
pkg/man/basic-probes.Rd
pkg/src/probe_acf.c
Log:
- add 'synth_loglik.c': synthetic likelihood computation a la Wood (2010)
- add 'correlation' option
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-10-18 17:41:01 UTC (rev 386)
+++ pkg/DESCRIPTION 2010-10-19 13:18:08 UTC (rev 387)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.34-3
-Date: 2010-10-18
+Date: 2010-10-19
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-10-18 17:41:01 UTC (rev 386)
+++ pkg/NAMESPACE 2010-10-19 13:18:08 UTC (rev 387)
@@ -18,6 +18,7 @@
probe_acf,
probe_ccf,
probe_nlar,
+ synth_loglik,
do_rprocess,
do_dprocess,
do_rmeasure,
Modified: pkg/R/basic-probes.R
===================================================================
--- pkg/R/basic-probes.R 2010-10-18 17:41:01 UTC (rev 386)
+++ pkg/R/basic-probes.R 2010-10-19 13:18:08 UTC (rev 387)
@@ -110,12 +110,13 @@
probe.acf <- function (var, lags, type = c("covariance", "correlation"), transform = identity) {
type <- match.arg(type)
+ corr <- type=="correlation"
transform <- match.fun(transform)
- corr <- type=="correlation"
if (corr && any(lags==0)) {
warning("useless zero lag discarded in ",sQuote("probe.acf"))
lags <- lags[lags!=0]
}
+ lags <- as.integer(lags)
function (y) .Call(
probe_acf,
x=transform(y[var,,drop=FALSE]),
@@ -124,15 +125,19 @@
)
}
-probe.ccf <- function (vars, lags, transform = identity) {
+probe.ccf <- function (vars, lags, type = c("covariance", "correlation"), transform = identity) {
+ type <- match.arg(type)
+ corr <- type=="correlation"
transform <- match.fun(transform)
if (length(vars)!=2)
stop(sQuote("vars")," must name two variables")
+ lags <- as.integer(lags)
function (y) .Call(
probe_ccf,
x=transform(y[vars[1],,drop=TRUE]),
y=transform(y[vars[2],,drop=TRUE]),
- lags=lags
+ lags=lags,
+ corr=corr
)
}
Modified: pkg/man/basic-probes.Rd
===================================================================
--- pkg/man/basic-probes.Rd 2010-10-18 17:41:01 UTC (rev 386)
+++ pkg/man/basic-probes.Rd 2010-10-19 13:18:08 UTC (rev 387)
@@ -24,7 +24,8 @@
probe.nlar(var, lags, powers, transform = identity)
probe.acf(var, lags, type = c("covariance", "correlation"),
transform = identity)
-probe.ccf(vars, lags, transform = identity)
+probe.ccf(vars, lags, type = c("covariance", "correlation"),
+ transform = identity)
probe.period(var, kernel.width, transform = identity)
probe.quantile(var, prob, transform = identity)
}
@@ -109,11 +110,14 @@
}
\item{\code{probe.acf}}{
returns a function that,
- if \code{type=="covariance"}, computes the autocovariance function of variable \code{var} at lags \code{lags};
- if \code{type=="correlation"}, computes the autocorrelation function of variable \code{var} at lags \code{lags}.
+ if \code{type=="covariance"}, computes the autocovariance of variable \code{var} at lags \code{lags};
+ if \code{type=="correlation"}, computes the autocorrelation of variable \code{var} at lags \code{lags}.
}
- \item{\code{probe.cov}, \code{probe.cor}}{
- 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.ccf}}{
+ returns a function that,
+ if \code{type=="covariance"}, computes the cross covariance of the two variables named in \code{vars} at lags \code{lags};
+ if \code{type=="correlation"}, computes the cross correlation.
+ }
\item{\code{probe.quantile}}{
returns a function that estimates the \code{prob}-th quantile of variable \code{var}.
}
Modified: pkg/src/probe_acf.c
===================================================================
--- pkg/src/probe_acf.c 2010-10-18 17:41:01 UTC (rev 386)
+++ pkg/src/probe_acf.c 2010-10-19 13:18:08 UTC (rev 387)
@@ -96,15 +96,15 @@
SEXP probe_acf (SEXP x, SEXP lags, SEXP corr) {
int nprotect = 0;
SEXP ans, ans_names;
- SEXP X, X_names, cov;
+ SEXP X, X_names;
int nlag, correlation, nvars, n;
int j, k, l;
- double *p, *p1;
+ double *p, *p1, *cov;
int *lag;
char tmp[BUFSIZ], *nm;
+ nlag = LENGTH(lags); // number of lags
PROTECT(lags = AS_INTEGER(lags)); nprotect++;
- nlag = LENGTH(lags); // number of lags
lag = INTEGER(lags);
correlation = *(INTEGER(AS_INTEGER(corr))); // correlation, or covariance?
@@ -120,9 +120,9 @@
if (correlation) {
l = 0;
- PROTECT(cov = NEW_NUMERIC(nvars)); nprotect++;
- pomp_acf_compute(REAL(cov),REAL(X),n,nvars,&l,1); // compute lag-0 covariance
- for (j = 0, p = REAL(ans), p1 = REAL(cov); j < nvars; j++, p1++)
+ cov = (double *) R_alloc(nvars,sizeof(double));
+ pomp_acf_compute(cov,REAL(X),n,nvars,&l,1); // compute lag-0 covariance
+ for (j = 0, p = REAL(ans), p1 = cov; j < nvars; j++, p1++)
for (k = 0; k < nlag; k++, p++)
*p /= *p1;
}
@@ -141,16 +141,18 @@
return ans;
}
-SEXP probe_ccf (SEXP x, SEXP y, SEXP lags) {
+SEXP probe_ccf (SEXP x, SEXP y, SEXP lags, SEXP corr) {
int nprotect = 0;
- SEXP ccf, ccf_names;
+ SEXP ans, ans_names;
SEXP X, Y;
- int nlag, n;
+ double cov[2], xx, *p;
+ int nlag, n, correlation;
int k;
char tmp[BUFSIZ], *nm;
nlag = LENGTH(lags);
PROTECT(lags = AS_INTEGER(lags)); nprotect++;
+ correlation = *(INTEGER(AS_INTEGER(corr))); // correlation, or covariance?
n = LENGTH(x); // n = # of observations
if (n != LENGTH(y))
@@ -159,18 +161,26 @@
PROTECT(X = duplicate(AS_NUMERIC(x))); nprotect++;
PROTECT(Y = duplicate(AS_NUMERIC(y))); nprotect++;
- PROTECT(ccf = NEW_NUMERIC(nlag)); nprotect++;
+ PROTECT(ans = NEW_NUMERIC(nlag)); nprotect++;
- pomp_ccf_compute(REAL(ccf),REAL(X),REAL(Y),n,INTEGER(lags),nlag);
+ pomp_ccf_compute(REAL(ans),REAL(X),REAL(Y),n,INTEGER(lags),nlag);
- PROTECT(ccf_names = NEW_STRING(nlag)); nprotect++;
+ if (correlation) {
+ k = 0;
+ pomp_acf_compute(&cov[0],REAL(X),n,1,&k,1); // compute lag-0 covariance of x
+ pomp_acf_compute(&cov[1],REAL(Y),n,1,&k,1); // compute lag-0 covariance of y
+ xx = sqrt(cov[0]*cov[1]);
+ for (k = 0, p = REAL(ans); k < nlag; k++, p++) *p /= xx; // convert to correlation
+ }
+
+ PROTECT(ans_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_STRING_ELT(ans_names,k,mkChar(tmp));
}
- SET_NAMES(ccf,ccf_names);
+ SET_NAMES(ans,ans_names);
UNPROTECT(nprotect);
- return ccf;
+ return ans;
}
Added: pkg/src/synth_lik.c
===================================================================
--- pkg/src/synth_lik.c (rev 0)
+++ pkg/src/synth_lik.c 2010-10-19 13:18:08 UTC (rev 387)
@@ -0,0 +1,117 @@
+// -*- mode: C++ -*-
+
+#include "pomp_internal.h"
+#include "pomp_mat.h"
+#include <string.h>
+
+// Campbell's robust variance-covariance estimation approach
+// described on p. 231 of Krzanowski 1988
+// with additional pre-conditioning for numerical stability
+// translated into C from R code by Simon N. Wood
+
+void robust_synth_loglik (double *y, int *dim, double *ydat, double *loglik) {
+ // 'ydat' is destroyed
+ // 'y' is preserved
+ int nrow = dim[0];
+ int ncol = dim[1];
+ double alpha = 2.0, beta = 1.25;
+ double x, xx, *yp, wbar;
+ double w[nrow], tau[ncol];
+ int pivot[nrow];
+ double *y1, *y2;
+ double d, d0, rss, half_log_det;
+ int i, j, k;
+
+ half_log_det = ncol*M_LN_SQRT_2PI;
+
+ y1 = (double *) R_alloc(nrow*ncol,sizeof(double));
+ y2 = (double *) R_alloc(nrow*ncol,sizeof(double));
+
+ // compute column means, center each column, compute preconditioner
+ memcpy(y1,y,nrow*ncol*sizeof(double));
+ for (yp = y1, j = 0; j < ncol; j++, yp += nrow) {
+ for (xx = 0, i = 0; i < nrow; i++) xx += yp[i];
+ xx /= nrow;
+ for (i = 0; i < nrow; i++) yp[i] -= xx; // center the column
+ for (xx = 0, i = 0; i < nrow; i++) xx += yp[i]*yp[i];
+ x = sqrt(xx/nrow); // column SD
+ for (i = 0; i < nrow; i++) yp[i] /= x; // precondition
+ }
+
+ // do first QR decomposition & backsolve
+
+ memcpy(y2,y1,nrow*ncol*sizeof(double));
+ pomp_qr(y2,nrow,ncol,pivot,tau); // Q*R = Y*inv(diag(d))
+ pomp_qr_x_inv_r(y2,nrow,ncol,y1,nrow); // y1 <- y1 %*% inv(R)
+
+ // create Campbell weight vector
+ d0 = sqrt(ncol)+alpha/sqrt(2.0);
+ for (wbar = 0, i = 0; i < nrow; i++) {
+ for (xx = 0, j = 0; j < ncol; j++) {
+ x = y1[i+nrow*j];
+ xx += x*x;
+ }
+ d = sqrt((nrow-1)*xx); // Mahalonobis distance of row i
+ if (d > d0) {
+ x = d-d0;
+ xx = exp(-0.5*x*x/beta)*d0/d;
+ } else {
+ xx = 1.0;
+ }
+ w[i] = xx;
+ wbar += xx*xx;
+ }
+ wbar = sqrt(wbar-1);
+
+ // compute weighted column means, center each column, precondition
+ memcpy(y1,y,nrow*ncol*sizeof(double));
+ for (yp = y1, j = 0; j < ncol; j++, yp += nrow) {
+ for (xx = 0, x = 0, i = 0; i < nrow; i++) {
+ xx += w[i]*yp[i];
+ x += w[i];
+ }
+ xx /= x; // column mean
+ for (i = 0; i < nrow; i++) yp[i] -= xx; // center the column
+ ydat[j] -= xx; // subtract mean from realized probe
+ for (xx = 0, i = 0; i < nrow; i++) {
+ xx += yp[i]*yp[i];
+ yp[i] /= wbar;
+ }
+ x = sqrt(xx/(nrow-1)); // column SD
+ for (i = 0; i < nrow; i++) yp[i] *= (w[i]/x); // precondition & weight
+ ydat[j] /= x;
+ half_log_det += log(x); // sum up logs(diag(D))
+ }
+
+ // do second QR decomposition & backsolve
+ pomp_qr(y1,nrow,ncol,pivot,tau); // Q*R = diag(w)*Y*inv(diag(d))
+ pomp_backsolve(y1,nrow,ncol,ydat,1,"Upper","Transpose","Non-unit"); // ydat <- ydat*inv(R)
+
+ // compute residual sum of squares and add up logs of diag(R)
+ xx = 0;
+ for (yp = y1, rss = 0, i = nrow+1, j = 0; j < ncol; j++, yp += i) { // yp marches along the diagonal of R
+ x = ydat[j];
+ rss += x*x;
+ half_log_det += log(fabs(*yp)); // log(diag(R))
+ xx += log(fabs(*yp));
+ }
+
+ *loglik = -0.5*rss-half_log_det;
+}
+
+
+SEXP synth_loglik (SEXP ysim, SEXP ydat) {
+ int nprotect = 0;
+ SEXP loglik, dim, y;
+
+ PROTECT(y = duplicate(AS_NUMERIC(ydat))); nprotect++;
+ PROTECT(dim = GET_DIM(ysim)); nprotect++;
+ PROTECT(ysim = AS_NUMERIC(ysim)); nprotect++;
+ PROTECT(loglik = NEW_NUMERIC(1)); nprotect++;
+
+ robust_synth_loglik(REAL(ysim),INTEGER(dim),REAL(y),REAL(loglik));
+
+ UNPROTECT(nprotect);
+ return loglik;
+}
+
Property changes on: pkg/src/synth_lik.c
___________________________________________________________________
Added: svn:eol-style
+ native
More information about the pomp-commits
mailing list