[Pomp-commits] r289 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 4 00:41:11 CEST 2010
Author: kingaa
Date: 2010-08-04 00:41:09 +0200 (Wed, 04 Aug 2010)
New Revision: 289
Added:
pkg/R/basic-probes.R
pkg/R/probe-match.R
pkg/R/probe.R
pkg/R/spect-match.R
pkg/R/spect.R
pkg/man/basic-probes.Rd
pkg/man/probe.Rd
pkg/man/probed-pomp-class.Rd
pkg/man/probed-pomp-methods.Rd
pkg/man/spect-pomp-class.Rd
pkg/man/spect.Rd
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
Log:
- add in the probe- and spectral-matching codes (thanks to Dan Reuman and Bruce Kendall)
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-08-03 20:42:25 UTC (rev 288)
+++ pkg/DESCRIPTION 2010-08-03 22:41:09 UTC (rev 289)
@@ -1,10 +1,10 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.31-2
-Date: 2010-07-19
+Version: 0.32-1
+Date: 2010-08-03
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing,
- Matthew J. Ferrari, Michael Lavine
+ Matthew J. Ferrari, Michael Lavine, Dan Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
URL: http://pomp.r-forge.r-project.org
Description: Inference methods for partially-observed Markov processes
@@ -21,3 +21,4 @@
mif-class.R particles-mif.R pfilter-mif.R mif.R mif-methods.R compare.mif.R
pmcmc.R pmcmc-methods.R compare.pmcmc.R
nlf-funcs.R nlf-guts.R nlf-lql.R nlf-objfun.R nlf.R
+ probe.R probe-match.R basic-probes.R spect.R spect-match.R
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-08-03 20:42:25 UTC (rev 288)
+++ pkg/NAMESPACE 2010-08-03 22:41:09 UTC (rev 289)
@@ -23,7 +23,7 @@
importFrom(deSolve,lsoda)
importFrom(subplex,subplex)
-exportClasses('pomp','mif','pmcmc')
+exportClasses('pomp','mif','pmcmc','probed.pomp','probe.matched.pomp','spect.pomp','spect.matched.pomp')
exportMethods(
'plot','show','print','coerce',
@@ -33,7 +33,8 @@
'particles','mif','continue','coef<-','states','trajectory',
'pred.mean','pred.var','filter.mean','conv.rec',
'bsmc',
- 'pmcmc','dprior'
+ 'pmcmc','dprior',
+ 'spect','probe'
)
export(
@@ -55,5 +56,15 @@
periodic.bspline.basis,
compare.mif,
nlf,
- traj.match
+ traj.match,
+ probe.mean,
+ probe.var,
+ probe.sd,
+ probe.period,
+ probe.quantile,
+ probe.acf,
+ probe.cov,
+ probe.cor,
+ probe.match,
+ spect.match
)
Added: pkg/R/basic-probes.R
===================================================================
--- pkg/R/basic-probes.R (rev 0)
+++ pkg/R/basic-probes.R 2010-08-03 22:41:09 UTC (rev 289)
@@ -0,0 +1,111 @@
+probe.mean <- function (var, trim = 0, transform = identity, na.rm = TRUE) {
+ function(y) mean(x=transform(y[var,]),trim=trim,na.rm=na.rm)
+}
+
+probe.var <- function (var, transform = identity, na.rm = TRUE) {
+ function(y) var(x=transform(y[var,]),na.rm=na.rm)
+}
+
+probe.sd <- function (var, transform = identity, na.rm = TRUE) {
+ function(y) sd(x=transform(y[var,]),na.rm=na.rm)
+}
+
+probe.period <- function (var, kernel.width, transform = identity) {
+ function (y) {
+ zz <- spec.pgram(
+ x=transform(y[var,]),
+ kernel=kernel("modified.daniell",m=kernel.width),
+ taper=0,
+ fast=FALSE,
+ pad=0,
+ detrend=FALSE,
+ plot=FALSE
+ )
+ 1/zz$freq[which.max(zz$spec)]
+ }
+}
+
+probe.quantile <- function (var, prob, transform = identity) {
+ function (y) quantile(transform(y[var,]),probs=prob)
+}
+
+probe.acf <- function (
+ var,
+ lag = 0,
+ type = c("correlation", "covariance", "partial"),
+ transform = identity,
+ ...
+ ) {
+ args <- list(...)
+ function (y) {
+ zz <- do.call(acf,c(list(x=transform(y[var,]),lag.max=lag,plot=FALSE),args))
+ if (type=="partial")
+ val <- zz$acf[lag]
+ else
+ val <- zz$acf[lag+1]
+ val
+ }
+}
+
+probe.cov <- function (
+ vars,
+ lag,
+ method = c("pearson", "kendall", "spearman"),
+ transform = identity
+ ) {
+ method <- match.arg(method)
+ lag <- as.integer(lag)
+ var1 <- vars[1]
+ if (length(vars)>1)
+ var2 <- vars[2]
+ else
+ var2 <- var1
+ function (y) {
+ if (lag>=0) {
+ val <- cov(
+ x=transform(y[var1,seq(from=1+lag,to=ncol(y),by=1)]),
+ y=transform(y[var2,seq(from=1,to=ncol(y)-lag,by=1)]),
+ method=method
+ )
+ } else {
+ val <- cov(
+ x=transform(y[var1,seq(from=1,to=ncol(y)+lag,by=1)]),
+ y=transform(y[var2,seq(from=-lag,to=ncol(y),by=1)]),
+ method=method
+ )
+ }
+ val
+ }
+}
+
+probe.cor <- function (
+ vars,
+ lag,
+ method = c("pearson", "kendall", "spearman"),
+ transform = identity
+ ) {
+ method <- match.arg(method)
+ lag <- as.integer(lag)
+ var1 <- vars[1]
+ if (length(vars)>1)
+ var2 <- vars[2]
+ else
+ var2 <- var1
+ function (y) {
+ if (lag>=0) {
+ val <- cor(
+ x=transform(y[var1,seq(from=1+lag,to=ncol(y),by=1)]),
+ y=transform(y[var2,seq(from=1,to=ncol(y)-lag,by=1)]),
+ method=method
+ )
+ } else {
+ val <- cor(
+ x=transform(y[var1,seq(from=1,to=ncol(y)+lag,by=1)]),
+ y=transform(y[var2,seq(from=-lag,to=ncol(y),by=1)]),
+ method=method
+ )
+ }
+ val
+ }
+}
+
Added: pkg/R/probe-match.R
===================================================================
--- pkg/R/probe-match.R (rev 0)
+++ pkg/R/probe-match.R 2010-08-03 22:41:09 UTC (rev 289)
@@ -0,0 +1,163 @@
+setClass(
+ "probe.matched.pomp",
+ contains="probed.pomp",
+ representation=representation(
+ weights="numeric",
+ fail.value="numeric",
+ evals="integer",
+ value="numeric",
+ convergence="integer",
+ msg="character"
+ )
+ )
+
+setMethod(
+ "summary",
+ "probe.matched.pomp",
+ function (object, ...) {
+ c(
+ summary(as(object,"probed.pomp")),
+ list(
+ weights=object at weights,
+ value=object at value,
+ eval=object at evals,
+ convergence=object at convergence
+ ),
+ if(length(object at msg)>0) list(msg=object at msg) else NULL
+ )
+ }
+ )
+
+probe.mismatch <- function (par, est, object, probes, params,
+ nsim = 1, seed = NULL,
+ weights, datval,
+ fail.value = NA) {
+ if (missing(par)) par <- numeric(0)
+ if (missing(est)) est <- integer(0)
+ if (missing(params)) params <- coef(object)
+
+ params[est] <- par
+
+ simval <- apply.probe.sim(object,probes=probes,params=params,nsim=nsim,seed=seed) # apply probes to model simulations
+
+ ## compute a measure of the discrepancies between simulations and data
+ sim.means <- colMeans(simval)
+ simval <- sweep(simval,2,sim.means)
+ discrep <- ((datval-sim.means)^2)/colMeans(simval^2)
+
+ if (!all(is.finite(discrep))) {
+ mismatch <- fail.value
+ } else {
+ mismatch <- sum(discrep*weights)/sum(weights)
+ }
+
+ mismatch
+}
+
+probe.match <- function(object, start, est = character(0),
+ probes, weights,
+ nsim, seed = NULL,
+ method = c("subplex","Nelder-Mead","SANN"),
+ verbose = getOption("verbose"),
+ eval.only = FALSE, fail.value = NA, ...) {
+
+ if (!is(object,"pomp"))
+ stop(sQuote("object")," must be of class ",sQuote("pomp"))
+
+ if (missing(start))
+ start <- coef(object)
+
+ if (!eval.only&&(length(est)<1))
+ stop("parameters to be estimated must be specified in ",sQuote("est"))
+ if (!is.character(est)|!all(est%in%names(start)))
+ stop(sQuote("est")," must refer to parameters named in ",sQuote("start"))
+ par.index <- which(names(start)%in%est)
+
+ if (missing(probes)) {
+ if (is(object,"probed.pomp"))
+ probes <- object at probes
+ else
+ stop(sQuote("probes")," must be supplied")
+ }
+ if (!is.list(probes)) probes <- list(probes)
+ if (!all(sapply(probes,is.function)))
+ stop(sQuote("probes")," must be a function or a list of functions")
+
+ if (missing(weights)) weights <- rep(1,length(probes))
+
+ method <- match.arg(method)
+
+ params <- start
+ guess <- params[par.index]
+
+ datval <- apply.probe.data(object,probes=probes) # apply probes to data
+
+ if (eval.only) {
+ val <- probe.mismatch(
+ par=guess,
+ est=par.index,
+ object=object,
+ probes=probes,
+ params=params,
+ nsim=nsim,
+ seed=seed,
+ weights=weights,
+ datval=datval,
+ fail.value=fail.value
+ )
+ conv <- 0
+ evals <- as.integer(c(1,0))
+ msg <- paste(sQuote("probe.mismatch"),"evaluated")
+ } else {
+ if (method == 'subplex') {
+ opt <- subplex::subplex(
+ par=guess,
+ fn=probe.mismatch,
+ est=par.index,
+ object=object,
+ probes=probes,
+ params=params,
+ nsim=nsim,
+ seed=seed,
+ weights=weights,
+ datval=datval,
+ fail.value=fail.value,
+ control=list(...)
+ )
+ } else {
+ opt <- optim(
+ par=guess,
+ fn=probe.mismatch,
+ est=par.index,
+ object=object,
+ probes=probes,
+ params=params,
+ nsim=nsim,
+ seed=seed,
+ weights=weights,
+ datval=datval,
+ fail.value=fail.value,
+ method=method,
+ control=list(...)
+ )
+ }
+ val <- opt$value
+ params[par.index] <- opt$par
+ conv <- opt$convergence
+ evals <- opt$counts
+ msg <- opt$message
+ }
+
+ coef(object,names(params)) <- unname(params)
+
+ new(
+ "probe.matched.pomp",
+ probe(object,probes=probes,nsim=nsim,seed=seed),
+ weights=weights,
+ fail.value=as.numeric(fail.value),
+ value=val,
+ convergence=as.integer(conv),
+ evals=as.integer(evals),
+ msg=as.character(msg)
+ )
+}
Added: pkg/R/probe.R
===================================================================
--- pkg/R/probe.R (rev 0)
+++ pkg/R/probe.R 2010-08-03 22:41:09 UTC (rev 289)
@@ -0,0 +1,221 @@
+setClass(
+ "probed.pomp",
+ contains="pomp",
+ representation(
+ probes="list",
+ seed="integer",
+ datvals="numeric",
+ simvals="array",
+ quantiles="numeric",
+ pvals="numeric"
+ )
+ )
+
+setMethod(
+ "summary",
+ "probed.pomp",
+ function (object, ...) {
+ list(
+ coef=coef(object),
+ nsim=nrow(object at simvals),
+ quantiles=object at quantiles,
+ pvals=object at pvals
+ )
+ }
+ )
+
+setGeneric("probe",function(object,probes,...)standardGeneric("probe"))
+
+apply.probe.data <- function (object, probes) {
+ val <- vector(mode="list",length=length(probes))
+ names(val) <- names(probes)
+ yy <- data.array(object)
+ for (p in seq_along(probes)) {
+ f <- probes[[p]]
+ if (length(formals(f))>1)
+ stop("each element of ",sQuote("probes")," must be a function of a single argument")
+ vv <- f(yy)
+ val[[p]] <- vv
+ }
+ do.call(c,val)
+}
+
+apply.probe.sim <- function (object, probes, params, nsim, seed) {
+ if (!is.numeric(nsim)||(length(nsim)>1)||(nsim<=0))
+ stop(sQuote("nsim")," must be a positive integer")
+ nsim <- as.integer(nsim)
+
+ y <- simulate(
+ object,
+ nsim=nsim,
+ seed=seed,
+ params=params,
+ times=time(object,t0=TRUE),
+ obs=TRUE
+ )[,,-1,drop=FALSE]
+ nmy <- dimnames(y)
+ dy <- dim(y)
+ yy <- array(dim=dy[c(1,3)],dimnames=nmy[c(1,3)])
+
+ val <- vector(mode="list",length=length(probes))
+ names(val) <- names(probes)
+ for (s in seq_len(nsim)) {
+ yy[,] <- y[,s,]
+ for (p in seq_along(probes)) {
+ f <- probes[[p]]
+ vv <- f(yy)
+ if (s==1) {
+ val[[p]] <- array(dim=c(nsim,length(vv)))
+ } else {
+ if (length(vv)!=ncol(val[[p]]))
+ if (length(vv)!=ncol(val))
+ stop("the size of the vector returned by a probe must not vary")
+ }
+ val[[p]][s,] <- vv
+ }
+ }
+ do.call(cbind,val)
+}
+
+setMethod(
+ "probe",
+ signature(object="pomp"),
+ function (object, probes, nsim = 1, seed = NULL, ...) {
+ if (!is.list(probes)) probes <- list(probes)
+ if (!all(sapply(probes,is.function)))
+ stop(sQuote("probes")," must be a function or a list of functions")
+ if (is.null(seed)) {
+ if (exists('.Random.seed',where=.GlobalEnv)) {
+ seed <- get(".Random.seed",pos=.GlobalEnv)
+ }
+ }
+
+ ## apply probes to data
+ datval <- apply.probe.data(object,probes=probes)
+ ## apply probes to model simulations
+ simval <- apply.probe.sim(
+ object,
+ probes=probes,
+ params=coef(object),
+ nsim=nsim,
+ seed=seed
+ )
+
+ pvals <- numeric(length(probes))
+ names(pvals) <- names(probes)
+ quants <- numeric(length(probes))
+ names(quants) <- names(probes)
+ for (k in seq_along(probes)) {
+ tails <- c(sum(simval[,k]>datval[k]),sum(simval[,k]<datval[k])+1)/(nsim+1)
+ pvals[k] <- min(c(2*tails,1))
+ quants[k] <- sum(simval[,k]<datval[k])/nsim
+ }
+
+ new(
+ "probed.pomp",
+ object,
+ probes=probes,
+ seed=as.integer(seed),
+ datvals=datval,
+ simvals=simval,
+ quantiles=quants,
+ pvals=pvals
+ )
+ }
+ )
+
+setMethod(
+ "probe",
+ signature(object="probed.pomp"),
+ function (object, probes, nsim, seed = NULL, ...) {
+ if (missing(probes)) probes <- object at probes
+ if (!is.list(probes)) probes <- list(probes)
+ if (!all(sapply(probes,is.function)))
+ stop(sQuote("probes")," must be a function or a list of functions")
+ if (is.null(seed)) {
+ if (exists('.Random.seed',where=.GlobalEnv)) {
+ seed <- get(".Random.seed",pos=.GlobalEnv)
+ }
+ }
+ if (missing(nsim)) nsim <- nrow(object at simvals)
+ probe(
+ as(object,"pomp"),
+ probes=probes,
+ nsim=nsim,
+ seed=seed,
+ ...
+ )
+ }
+ )
+
+setMethod(
+ "plot",
+ "probed.pomp",
+ function (x, y, ...) {
+
+ ##function for plotting diagonal panels
+ diag.panel.hist <- function(x, ...) {
+ ##plot a histogram for the simulations
+ usr <- par("usr")
+ on.exit(par(usr))
+ par(usr=c(usr[1:2],0,1.5))
+ h <- hist(x[-1],plot=FALSE)
+ breaks <- h$breaks
+ nB <- length(breaks)
+ y <- h$counts
+ y <- y/max(y)
+ rect(breaks[-nB],0,breaks[-1],y,...)
+ ##plot the data point
+ lines(c(x[1],x[1]),c(0,max(h$counts)),col="red")
+ }
+
+ ##function for plotting above-diagonal panels
+ above.diag.panel <- function (x, y, ...) {
+ ##plot the simulations
+ points(x[-1],y[-1],...)
+ ##plot the data
+ mMx <- c(min(x),max(x))
+ mMy <- c(min(y),max(y))
+ lines(c(x[1],x[1]),mMy,col="red")
+ lines(mMx,c(y[1],y[1]),col="red")
+ }
+
+ ##function for plotting below-diagonal panels
+ below.diag.panel <- function (x, y, ...) {
+ mMx <- c(min(x),max(x))
+ mMy <- c(min(y),max(y))
+ x <- x[-1]
+ y <- y[-1]
+ correls <- round(cor(x,y),3)
+ text(mean(mMx),mean(mMy),correls,cex=1)
+ }
+
+ ##prepare the arguments for these functions
+ nprobes <- length(x at datvals)
+ nsim <- nrow(x at simvals)
+ datsimvals <- array(dim=c(nsim+1,nprobes))
+ datsimvals[1,] <- x at datvals
+ datsimvals[-1,] <- x at simvals
+
+ labels <- paste("pb",seq_len(nprobes))
+ if (!is.null(names(x at datvals)))
+ labels <- ifelse(names(x at datvals)=="",labels,names(x at datvals))
+ lab.plus <- paste(labels,paste("p=",round(x at pvals,3),sep=""),sep="\n")
+ ##now make the plot
+
+ if (nprobes>1) {
+ pairs(
+ datsimvals,
+ diag.panel=diag.panel.hist,
+ lower.panel=below.diag.panel,
+ upper.panel=above.diag.panel,
+ labels=lab.plus,
+ cex.labels=1
+ )
+ } else {
+ plot(datsimvals,datsimvals,type="n",xlab="",ylab="",yaxt="n",main=lab.plus)
+ diag.panel.hist(datsimvals)
+ }
+ }
+ )
+
Added: pkg/R/spect-match.R
===================================================================
--- pkg/R/spect-match.R (rev 0)
+++ pkg/R/spect-match.R 2010-08-03 22:41:09 UTC (rev 289)
@@ -0,0 +1,244 @@
+setClass(
+ "spect.matched.pomp",
+ contains="spect.pomp",
+ representation=representation(
+ est="character",
+ fail.value="numeric",
+ evals="integer",
+ value="numeric",
+ weights="numeric",
+ convergence="integer",
+ msg="character"
+ )
+ )
+
+setMethod(
+ "summary",
+ "spect.matched.pomp",
+ function (object, ...) {
+ c(
+ summary(as(object,"spect.pomp")),
+ list(
+ est=object at est,
+ value=object at value,
+ eval=object at evals,
+ convergence=object at convergence
+ ),
+ if(length(object at msg)>0) list(msg=object at msg) else NULL
+ )
+ }
+ )
+
+spect.mismatch <- function (par, est, object, params,
+ vars, ker, nsim, seed,
+ transform, detrend, weights,
+ data.spec, fail.value) {
+ if (missing(par)) par <- numeric(0)
+ if (missing(est)) est <- integer(0)
+ if (missing(params)) params <- coef(object)
+
+ params[est] <- par
+
+ ## vector of frequencies and estimated power spectum of data
+ freq <- data.spec$freq
+ datval <- data.spec$spec
+
+ if (missing(weights)) weights <- 1
+
+ ## estimate power spectra of simulations
+ simvals <- compute.spect.sim(
+ object,
+ vars=vars,
+ params=params,
+ nsim=nsim,
+ seed=seed,
+ transform=transform,
+ detrend=detrend,
+ ker=ker
+ )
+
+ ## compute a measure of the discrepancies between simulations and data
+ discrep <- array(dim=c(length(freq),length(vars)))
+ sim.means <- colMeans(simvals)
+ for (j in seq_along(freq)) {
+ for (k in seq_along(vars)) {
+ discrep[j,k] <- ((datval[j,k]-sim.means[j,k])^2)/mean((simvals[j,,k]-sim.means[j,k])^2)
+ }
+ discrep[j,] <- weights[j]*discrep[j,]
+ }
+
+ if (!all(is.finite(discrep))) {
+ mismatch <- fail.value
+ } else {
+ mismatch <- sum(discrep)
+ }
+
+ mismatch
+}
+
+spect.match <- function(object, start, est = character(0),
+ vars, nsim, seed = NULL,
+ kernel.width, transform = identity,
+ detrend = c("none","mean","linear","quadratic"),
+ weights,
+ method = c("subplex","Nelder-Mead","SANN"),
+ verbose = getOption("verbose"),
+ eval.only = FALSE, fail.value = NA, ...) {
+
+ if (!is(object,"pomp"))
+ stop(sQuote("object")," must be of class ",sQuote("pomp"))
+
+ if (missing(vars))
+ vars <- rownames(object at data)
+
+ if (missing(kernel.width)) {
+ if (is(object,"spect.pomp")) {
+ kernel.width <- object at kernel.width
+ } else {
+ stop(sQuote("kernel.width")," must be specified")
+ }
+ }
+ ker <- reuman.kernel(kernel.width)
+
+ if (missing(transform)) {
+ if (is(object,"spect.pomp")) {
+ transform <- object at transform
+ } else {
+ transform <- identity
+ }
+ }
+
+ if (missing(nsim)||(nsim<1))
+ stop(sQuote("nsim")," must be specified as a positive integer")
+
+ if (missing(detrend)) {
+ if (is(object,"spect.pomp")) {
+ detrend <- object at detrend
+ } else {
+ detrend <- "none"
+ }
+ }
+ detrend <- match.arg(detrend)
+
+ method <- match.arg(method)
+
+ ds <- compute.spect.data(
+ object,
+ vars=vars,
+ transform=transform,
+ detrend=detrend,
+ ker=ker
+ )
+
+ if (missing(weights)) weights <- 1
+ if (is.numeric(weights)) {
+ if ((length(weights)!=1)&&(length(weights)!=length(ds$freq)))
+ stop("if ",sQuote("weights")," is provided as a vector, it must have length ",length(ds$freq))
+ } else if (is.function(weights)) {
+ weights <- sapply(ds$freq,weights)
+ } else {
+ stop(sQuote("weights")," must be specified as a vector or as a function")
+ }
+ if (any((!is.finite(weights))|(weights<0)))
+ stop(sQuote("weights")," should be nonnegative and finite")
+ weights <- weights/mean(weights)
+
+ if (missing(start))
+ start <- coef(object)
+
+ if (!eval.only&&(length(est)<1))
+ stop("parameters to be estimated must be specified in ",sQuote("est"))
+ if (!is.character(est)|!all(est%in%names(start)))
+ stop(sQuote("est")," must refer to parameters named in ",sQuote("start"))
+ par.index <- which(names(start)%in%est)
+
+ params <- start
+ guess <- params[par.index]
+
+ if (eval.only) {
+ val <- spect.mismatch(
+ par=guess,
+ est=par.index,
+ object=object,
+ params=params,
+ vars=vars,
+ ker=ker,
+ nsim=nsim,
+ seed=seed,
+ transform=transform,
+ detrend=detrend,
+ weights=weights,
+ data.spec=ds,
+ fail.value=fail.value
+ )
+ conv <- 0
+ evals <- as.integer(1)
+ msg <- paste(sQuote("spec.mismatch"),"evaluated")
+ } else {
+ if (method == 'subplex') {
+ opt <- subplex::subplex(
+ par=guess,
+ fn=spect.mismatch,
+ est=par.index,
+ object=object,
+ params=params,
+ vars=vars,
+ ker=ker,
+ nsim=nsim,
+ seed=seed,
+ transform=transform,
+ detrend=detrend,
+ weights=weights,
+ data.spec=ds,
+ fail.value=fail.value,
+ control=list(...)
+ )
+ } else {
+ opt <- optim(
+ par=guess,
+ fn=spect.mismatch,
+ est=par.index,
+ object=object,
+ params=params,
+ vars=vars,
+ ker=ker,
+ nsim=nsim,
+ seed=seed,
+ transform=transform,
+ detrend=detrend,
+ weights=weights,
+ data.spec=ds,
+ fail.value=fail.value,
+ method=method,
+ control=list(...)
+ )
+ }
+ val <- opt$value
+ params[par.index] <- opt$par
+ conv <- opt$convergence
+ evals <- opt$counts
+ msg <- opt$message
+ }
+
+ coef(object,names(params)) <- unname(params)
+
+ new(
+ "spect.matched.pomp",
+ spect(
+ object,
+ vars=vars,
+ kernel.width=kernel.width,
+ nsim=nsim,
+ seed=seed,
+ transform=transform,
+ detrend=detrend
+ ),
+ est=names(start)[par.index],
+ fail.value=as.numeric(fail.value),
+ value=val,
+ weights=weights,
+ convergence=as.integer(conv),
+ evals=as.integer(evals),
+ msg=as.character(msg)
+ )
+}
Added: pkg/R/spect.R
===================================================================
--- pkg/R/spect.R (rev 0)
+++ pkg/R/spect.R 2010-08-03 22:41:09 UTC (rev 289)
@@ -0,0 +1,362 @@
+# Authors:
+# Cai GoGwilt, MIT
+# Daniel Reuman, Imperial College London
+# Aaron A. King, U of Michigan
+
+setClass(
+ "spect.pomp",
+ contains="pomp",
+ representation=representation(
+ seed="integer",
+ kernel.width="numeric",
+ transform="function",
+ freq="numeric",
+ datspec="array",
+ simspec="array",
+ pvals="numeric",
+ detrend="character"
+ )
+ )
+
+setMethod(
+ "summary",
+ "spect.pomp",
+ function (object, ...) {
+ list(
+ coef=coef(object),
+ nsim=nrow(object at simspec),
+ pvals=object at pvals
+ )
+ }
+ )
+
+## detrends in one of several ways, according to type.
+## tseries is a numeric vector,
+pomp.detrend <- function (tseries, type)
+ switch(
+ type,
+ mean=tseries-mean(tseries),
+ linear={
+ x <- seq_along(tseries)
+ lm(tseries~x)$residuals
+ },
+ quadratic={
+ x <- seq_along(tseries)
+ lm(tseries~x+I(x^2))$residuals
+ },
+ tseries
+ )
+
+## The default smoothing kernel for the R spec.pgram function is weird.
+## This function creates a better one.
+reuman.kernel <- function (kernel.width) {
+ ker <- kernel("modified.daniell",m=kernel.width)
+ x <- seq.int(from=0,to=kernel.width,by=1)/kernel.width
+ ker[[1]] <- (15/(16*2*pi))*((x-1)^2)*((x+1)^2)
+ ker[[1]] <- ker[[1]]/(2*sum(ker[[1]][-1])+ker[[1]][1])
+ attr(ker,"name") <- NULL
+ ker
+}
+
+compute.spect.data <- function (object, vars, transform, detrend, ker) {
+ dat <- data.array(object,vars)
+ if (any(is.na(dat)))
+ stop(sQuote("spect")," is incompatible with NAs in the data")
+ dt <- diff(time(object,t0=FALSE))
+ dt.tol <- 0.025
+ if (max(dt)-min(dt)>dt.tol*mean(dt))
+ stop(sQuote("spect")," assumes evenly spaced times")
+ for (j in seq_along(vars)) {
+ sp <- spec.pgram(
+ pomp.detrend(
+ transform(dat[j,]),
+ type=detrend
+ ),
+ spans=ker,
+ taper=0,
+ pad=0,
+ fast=FALSE,
+ detrend=FALSE,
+ plot=FALSE
+ )
+ if (j==1) {
+ freq <- sp$freq
+ datspec <- array(
+ dim=c(length(freq),nrow(dat)),
+ dimnames=list(NULL,vars)
+ )
+ }
+ datspec[,j] <- log10(sp$spec)
+ }
+ list(freq=freq,spec=datspec)
+}
+
+compute.spect.sim <- function (object, vars, params, nsim, seed, transform, detrend, ker) {
+ sims <- try(
+ simulate(
+ object,
+ nsim=nsim,
+ seed=seed,
+ params=params,
+ times=time(object,t0=TRUE),
+ obs=TRUE
+ )[,,-1,drop=FALSE],
+ silent=FALSE
+ )
+ if (inherits(sims,"try-error"))
+ stop(sQuote("spect")," error: cannot simulate")
+ sims <- sims[vars,,,drop=FALSE]
+ if (any(is.na(sims)))
+ stop("NA in simulated data series")
+ nobs <- length(vars)
+ for (j in seq_len(nobs)) {
+ for (k in seq_len(nsim)) {
+ sp <- spec.pgram(
+ pomp.detrend(
+ transform(sims[j,k,]),
+ type=detrend
+ ),
+ spans=ker,
+ taper=0,
+ pad=0,
+ fast=FALSE,
+ detrend=FALSE,
+ plot=FALSE
+ )
+ if ((j==1)&&(k==1)) {
+ simspec <- array(
+ dim=c(nsim,length(sp$freq),nobs),
+ dimnames=list(NULL,NULL,vars)
+ )
+ }
+ simspec[k,,j] <- log10(sp$spec)
+ }
+ }
+ simspec
+}
+
+setGeneric("spect",function(object,...)standardGeneric("spect"))
+
+setMethod(
+ "spect",
+ signature(object="pomp"),
+ function (object, vars, kernel.width, nsim, seed = NULL, transform = identity,
+ detrend = c("none","mean","linear","quadratic"),
+ ...) {
+
+ if (missing(vars))
+ vars <- rownames(object at data)
+
+ if (missing(kernel.width))
+ stop(sQuote("kernel.width")," must be specified")
+ if (missing(nsim)||(nsim<1))
+ stop(sQuote("nsim")," must be specified as a positive integer")
+
+ detrend <- match.arg(detrend)
+ ker <- reuman.kernel(kernel.width)
+
+ if (is.null(seed)) {
+ if (exists('.Random.seed',where=.GlobalEnv)) {
+ seed <- get(".Random.seed",pos=.GlobalEnv)
+ }
+ }
+
+ ds <- compute.spect.data(
+ object,
+ vars=vars,
+ transform=transform,
+ detrend=detrend,
+ ker=ker
+ )
+ freq <- ds$freq
+ datspec <- ds$spec
+ simspec <- compute.spect.sim(
+ object,
+ vars=vars,
+ params=coef(object),
+ nsim=nsim,
+ seed=seed,
+ transform=transform,
+ detrend=detrend,
+ ker=ker
+ )
+
+ pvals <- numeric(length(vars)+1)
+ names(pvals) <- c(vars,"all")
+ mean.simspec <- colMeans(simspec) # mean spectrum of simulations
+ totdatdist <- 0
+ totsimdist <- 0
+ for (j in seq_along(vars)) {
+ ## L-2 distance between data and mean simulated spectrum
+ datdist <- sum((datspec[,j]-mean.simspec[,j])^2)
+ ## L-2 distance betw. each sim. and mean simulated spectrum
+ simdist <- sapply(
+ seq_len(nsim),
+ function(k)sum((simspec[k,,j]-mean.simspec[,j])^2)
+ )
+ pvals[j] <- (nsim+1-sum(simdist<datdist))/(nsim+1)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 289
More information about the pomp-commits
mailing list