[Pomp-commits] r214 - in pkg: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Apr 29 21:54:45 CEST 2010
Author: kingaa
Date: 2010-04-29 21:54:45 +0200 (Thu, 29 Apr 2010)
New Revision: 214
Modified:
pkg/R/nlf-funcs.R
pkg/R/nlf-lql.R
pkg/R/nlf-objfun.R
pkg/R/nlf.R
pkg/man/dmeasure-pomp.Rd
pkg/man/dprocess-pomp.Rd
pkg/man/init.state-pomp.Rd
pkg/man/mif-class.Rd
pkg/man/nlf.Rd
pkg/man/particles-mif.Rd
pkg/man/pomp-class.Rd
pkg/man/rmeasure-pomp.Rd
pkg/man/rprocess-pomp.Rd
pkg/man/skeleton-pomp.Rd
Log:
- nlf now takes an optional argument 'transform' which transforms the data before fitting
- the low-level interface documentation now has keyword 'internal' which keeps it out of the manual
Modified: pkg/R/nlf-funcs.R
===================================================================
--- pkg/R/nlf-funcs.R 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/R/nlf-funcs.R 2010-04-29 19:54:45 UTC (rev 214)
@@ -65,7 +65,7 @@
}
}
-Newey.West <-function(x,y,maxlag) {
+Newey.West <- function(x,y,maxlag) {
w <- 1-(1:maxlag)/(maxlag+1)
out <- mean(x*y,na.rm=T)
for(i in 1:maxlag) {
@@ -75,16 +75,16 @@
}
-make.tensorbasis.NLF=function(A,B) {
- if(nrow(A)!=nrow(B)) stop("Incompatible matrices in make.tensorbasis");
- ncol.A=ncol(A)
- ncol.B=ncol(B)
- Tmat=matrix(0,nrow(A),ncol.A*ncol.B)
- for(i in 1:ncol.A) {
- start=(i-1)*ncol.B
- for(j in 1:ncol.B) {
- Tmat[,start+j]=A[,i]*B[,j]
- }
- }
- return(Tmat)
+make.tensorbasis.NLF <- function(A,B) {
+ if(nrow(A)!=nrow(B)) stop("Incompatible matrices in make.tensorbasis")
+ ncol.A <- ncol(A)
+ ncol.B <- ncol(B)
+ Tmat <- matrix(0,nrow(A),ncol.A*ncol.B)
+ for(i in 1:ncol.A) {
+ start=(i-1)*ncol.B
+ for(j in 1:ncol.B) {
+ Tmat[,start+j] <- A[,i]*B[,j]
+ }
+ }
+ Tmat
}
Modified: pkg/R/nlf-lql.R
===================================================================
--- pkg/R/nlf-lql.R 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/R/nlf-lql.R 2010-04-29 19:54:45 UTC (rev 214)
@@ -1,5 +1,6 @@
NLF.LQL <- function (params.fitted, object, params, par.index,
- times, lags, period, tensor, seed, nrbf = 4, verbose = FALSE,
+ times, lags, period, tensor, seed = NULL, transform = function(x)x,
+ nrbf = 4, verbose = FALSE,
bootstrap = FALSE, bootsamp = NULL) {
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -12,7 +13,6 @@
params[par.index] <- params.fitted
params <- as.matrix(params)
- xstart <- init.state(object,params)
## Need to extract number of state variables (nvar) from pomp object
## Need to include simulation times in problem specification
@@ -21,37 +21,20 @@
## Version 0.2, May 2008, Stephen P. Ellner
data.ts <- data.array(object)
-
- ## set the random seed (be very careful about this)
- if (!is.null(seed)) {
- if (!exists('.Random.seed',envir=.GlobalEnv)) runif(1)
- save.seed <- get('.Random.seed',envir=.GlobalEnv)
- set.seed(seed)
- }
-
- x <- try(
- rprocess(object,xstart=xstart,times=times,params=params), # simulate the process model
- silent=FALSE
- )
- if (inherits(x,'try-error'))
- stop("'NLF.LQL' reports: error in user 'rprocess'")
-
- if (!all(is.finite(x))) return(FAILED)
-
+
y <- try(
- rmeasure(object,x=x[,,-1,drop=F],times=times[-1],params=params),
+ simulate(object,times=times,params=params,seed=seed,obs=TRUE,states=FALSE),
silent=FALSE
)
- if (inherits(y,'try-error'))
- stop("'NLF.LQL' reports: error in user 'rmeasure'")
+ if (inherits(y,"try-error"))
+ stop(sQuote("NLF.LQL")," reports: error in simulation")
+ ## Test whether the model time series is valid
+ if (!all(is.finite(y))) return(FAILED)
- if (!is.null(seed)) assign('.Random.seed',save.seed,envir=.GlobalEnv) # restore the RNG state
-
- ## Test whether the model time series is valid
model.ts <- array(dim=c(nrow(data.ts),length(times)-1),dimnames=list(rownames(data.ts),NULL))
- if (!all(is.finite(x))) return(FAILED)
-
- model.ts[,] <- y[,1,]
+ model.ts[,] <- apply(y[,1,-1,drop=FALSE],c(2,3),transform)
+ data.ts[,] <- apply(data.ts,2,transform)
+
LQL <- try(
NLF.guts(
data.mat=data.ts,
@@ -59,17 +42,18 @@
model.mat=model.ts,
model.times=times[-1],
lags=lags,
- period=period, tensor=tensor,
+ period=period,
+ tensor=tensor,
nrbf=nrbf,
- verbose=F,
+ verbose=FALSE,
bootstrap,
bootsamp,
- plotfit=F
+ plotfit=FALSE
),
silent=FALSE
)
if (inherits(LQL,"try-error"))
- stop("'NLF.LQL' reports: error in 'NLF.guts'")
+ stop(sQuote("NLF.LQL")," reports: error in ",sQuote("NLF.guts"))
LQL
}
Modified: pkg/R/nlf-objfun.R
===================================================================
--- pkg/R/nlf-objfun.R 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/R/nlf-objfun.R 2010-04-29 19:54:45 UTC (rev 214)
@@ -1,21 +1,23 @@
nlf.objfun <- function (params.fitted, object, params, par.index,
- times, lags, period, tensor, seed, nrbf = 4,
- verbose = FALSE, bootstrap = FALSE, bootsamp = NULL)
+ times, lags, period, tensor, seed, transform = function(x)x,
+ nrbf = 4, verbose = FALSE, bootstrap = FALSE, bootsamp = NULL)
{
-sum(
NLF.LQL(
- params.fitted,
- object,
- params,
- par.index,
- times,
- lags,
- period, tensor,
- seed,
- nrbf,
- verbose,
- bootstrap,
- bootsamp
+ params.fitted=params.fitted,
+ object=object,
+ params=params,
+ par.index=par.index,
+ times=times,
+ lags=lags,
+ period=period,
+ tensor=tensor,
+ seed=seed,
+ transform=transform,
+ nrbf=nrbf,
+ verbose=verbose,
+ bootstrap=bootstrap,
+ bootsamp=bootsamp
),
na.rm=T
)
Modified: pkg/R/nlf.R
===================================================================
--- pkg/R/nlf.R 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/R/nlf.R 2010-04-29 19:54:45 UTC (rev 214)
@@ -1,5 +1,6 @@
nlf <- function(object, start, est, lags, period = NA, tensor = FALSE, nconverge = 1000, nasymp = 1000,
- seed = 1066, nrbf = 4, method = "subplex", skip.se = FALSE, verbose = FALSE, gr = NULL,
+ seed = 1066, transform = function(x)x,
+ nrbf = 4, method = "subplex", skip.se = FALSE, verbose = FALSE, gr = NULL,
bootstrap = FALSE, bootsamp = NULL, lql.frac = 0.1, se.par.frac = 0.1, eval.only = FALSE, ...) {
## Fit a POMP object using NLF
## v. 0.1, 3 Dec. 2007
@@ -19,6 +20,9 @@
if (!is(object,'pomp'))
stop("'object' must be a 'pomp' object")
+ if (!is.function(transform))
+ stop(sQuote("transform")," must be a function!")
+
if (eval.only) est <- 1
if (is.character(est)) {
@@ -60,29 +64,23 @@
if (eval.only) {
- opt <- optim(
- par=guess,
- fn=nlf.objfun,
- gr=gr,
- method="SANN",
- object=object,
- params=params,
- par.index=par.index,
- times=times,
- lags=lags,
- period=period,
- tensor=tensor,
- seed=seed,
- nrbf=nrbf,
- hessian=FALSE,
- verbose=verbose,
- bootstrap=bootstrap,
- bootsamp=bootsamp,
- control=list(
- maxit=0
- )
- )
- return(-opt$value)
+ val <- nlf.objfun(
+ params.fitted=guess,
+ object=object,
+ params=params,
+ par.index=par.index,
+ times=times,
+ lags=lags,
+ period=period,
+ tensor=tensor,
+ seed=seed,
+ transform=transform,
+ nrbf=nrbf,
+ verbose=verbose,
+ bootstrap=bootstrap,
+ bootsamp=bootsamp
+ )
+ return(-val)
}
if (method == 'subplex') {
@@ -97,6 +95,7 @@
period=period,
tensor=tensor,
seed=seed,
+ transform=transform,
nrbf=nrbf,
verbose=verbose,
bootstrap=bootstrap,
@@ -117,6 +116,7 @@
period=period,
tensor=tensor,
seed=seed,
+ transform=transform,
nrbf=nrbf,
verbose=verbose,
bootstrap=bootstrap,
@@ -136,7 +136,8 @@
Jhat <- matrix(0,nfitted,nfitted)
Ihat <- Jhat
f0 <- NLF.LQL(fitted,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed,
+ transform=transform, nrbf=4,
verbose=FALSE)
F0 <- mean(f0,na.rm=T)
@@ -154,22 +155,22 @@
guess <- fitted
guess[i] <- fitted[i]-sqrt(2)*h*abs(fitted[i])
Fvals[1] <- mean(NLF.LQL(guess,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
- verbose=FALSE),na.rm=T)
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed, transform=transform,
+ nrbf=4, verbose=FALSE),na.rm=T)
guess <- fitted
guess[i] <- fitted[i]-h*abs(fitted[i])
Fvals[2] <- mean(NLF.LQL(guess,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T);
guess <- fitted
guess[i] <- fitted[i]+h*abs(fitted[i])
Fvals[4] <- mean(NLF.LQL(guess,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
guess <- fitted
guess[i] <- fitted[i]+sqrt(2)*h*abs(fitted[i])
Fvals[5] <- mean(NLF.LQL(guess,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T);
FAILED = - 999999
Fvals[Fvals < FAILED+10] <- NA
@@ -185,12 +186,12 @@
guess.up <- fitted
guess.up[i] <- guess.up[i]+eps[i]
f.up <- NLF.LQL(guess.up,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed, transform=transform, nrbf=4,
verbose=FALSE)
F.up <- mean(f.up,na.rm=T)
f.up2 <- NLF.LQL(guess.up,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed, transform=transform, nrbf=4,
verbose=FALSE)
if (verbose) cat("Fitted param ", i, F.up, mean(f.up2,na.rm=T)," up in ",sQuote("nlf"),"\n")
@@ -198,7 +199,7 @@
guess.down <- fitted
guess.down[i] <- guess.down[i]-eps[i]
f.down <- NLF.LQL(guess.down,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed, transform=transform, nrbf=4,
verbose=FALSE)
F.down <- mean(f.down,na.rm=T)
@@ -215,28 +216,28 @@
guess.uu[i] <- guess.uu[i]+eps[i]
guess.uu[j] <- guess.uu[j]+eps[j]
F.uu <- mean(NLF.LQL(guess.uu,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
guess.ud <- fitted
guess.ud[i] <- guess.ud[i]+eps[i]
guess.ud[j] <- guess.ud[j]-eps[j]
F.ud <- mean(NLF.LQL(guess.ud,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
guess.du <- fitted
guess.du[i] <- guess.du[i]-eps[i]
guess.du[j] <- guess.du[j]+eps[j]
F.du <- mean(NLF.LQL(guess.du,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
guess.dd <- fitted
guess.dd[i] <- guess.dd[i]-eps[i]
guess.dd[j] <- guess.dd[j]-eps[j]
F.dd <- mean(NLF.LQL(guess.dd,object=object, params=params, par.index=par.index,
- times=times, lags=lags, period=period, tensor=tensor, seed=seed, nrbf=4,
+ times=times, lags=lags, period=period, tensor=tensor, seed=seed, transform=transform, nrbf=4,
verbose=FALSE),na.rm=T)
dij <- (F.uu+F.dd)-(F.ud+F.du)
Modified: pkg/man/dmeasure-pomp.Rd
===================================================================
--- pkg/man/dmeasure-pomp.Rd 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/man/dmeasure-pomp.Rd 2010-04-29 19:54:45 UTC (rev 214)
@@ -3,6 +3,7 @@
\alias{dmeasure}
\alias{dmeasure,pomp-method}
\alias{dmeasure-pomp}
+\keyword{internal}
\title{Evaluate the probability density of observations given underlying states in a partially-observed Markov process}
\description{
The method \code{dmeasure} evaluates the probability density of a set of measurements given the state of the system.
Modified: pkg/man/dprocess-pomp.Rd
===================================================================
--- pkg/man/dprocess-pomp.Rd 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/man/dprocess-pomp.Rd 2010-04-29 19:54:45 UTC (rev 214)
@@ -3,6 +3,7 @@
\alias{dprocess}
\alias{dprocess,pomp-method}
\alias{dprocess-pomp}
+\keyword{internal}
\title{
Evaluate the probability density of state transitions in a Markov process
}
Modified: pkg/man/init.state-pomp.Rd
===================================================================
--- pkg/man/init.state-pomp.Rd 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/man/init.state-pomp.Rd 2010-04-29 19:54:45 UTC (rev 214)
@@ -3,6 +3,7 @@
\alias{init.state}
\alias{init.state,pomp-method}
\alias{init.state-pomp}
+\keyword{internal}
\title{Return a matrix of initial conditions given a vector of parameters and an initial time.}
\description{
The method \code{init.state} returns a vector of initial conditions for the state process when given a vector of parameters \code{params} and an initial time \code{t0}.
Modified: pkg/man/mif-class.Rd
===================================================================
--- pkg/man/mif-class.Rd 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/man/mif-class.Rd 2010-04-29 19:54:45 UTC (rev 214)
@@ -1,6 +1,7 @@
\name{mif-class}
\docType{class}
\alias{mif-class}
+\keyword{internal}
\title{The "mif" class}
\description{
The \code{mif} class holds a fitted model and is created by a call to \code{\link{mif}}.
Modified: pkg/man/nlf.Rd
===================================================================
--- pkg/man/nlf.Rd 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/man/nlf.Rd 2010-04-29 19:54:45 UTC (rev 214)
@@ -7,8 +7,8 @@
}
\usage{
nlf(object, start, est, lags, period = NA, tensor = FALSE,
- nconverge=1000, nasymp=1000, seed = 1066, nrbf = 4,
- method = "subplex", skip.se = FALSE, verbose = FALSE, gr = NULL,
+ nconverge=1000, nasymp=1000, seed = 1066, transform = function (x) x,
+ nrbf = 4, method = "subplex", skip.se = FALSE, verbose = FALSE, gr = NULL,
bootstrap=FALSE, bootsamp = NULL,
lql.frac = 0.1, se.par.frac = 0.1, eval.only = FALSE, \dots)
}
@@ -46,6 +46,14 @@
When fitting, it is usually best to always run the simulations with the same sequence of random numbers, which is accomplished by setting \code{seed} to an integer.
If you want a truly random simulation, set \code{seed=NULL}.
}
+ \item{transform}{
+ optional function.
+ If specified, forecasting is performed using data and model simulations transformed by this function.
+ By default, \code{transform} is the identity function.
+ The main purpose of \code{transform} is to achieve approximately multivariate normal forecasting errors.
+ If data are univariate, \code{transform} should take a scalar and return a scalar.
+ If data are multivariate, \code{transform} should assume a vector input and return a vector of the same length.
+ }
\item{nrbf}{A scalar specifying the number of radial basis functions to be used at each lag.}
\item{method}{
Optimization method.
Modified: pkg/man/particles-mif.Rd
===================================================================
--- pkg/man/particles-mif.Rd 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/man/particles-mif.Rd 2010-04-29 19:54:45 UTC (rev 214)
@@ -3,6 +3,7 @@
\alias{particles}
\alias{particles,mif-method}
\alias{particles-mif}
+\keyword{internal}
\title{Generate particles from the user-specified distribution.}
\description{
Generate particles from the user-specified distribution.
Modified: pkg/man/pomp-class.Rd
===================================================================
--- pkg/man/pomp-class.Rd 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/man/pomp-class.Rd 2010-04-29 19:54:45 UTC (rev 214)
@@ -1,7 +1,8 @@
\name{pomp-class}
\docType{class}
\alias{pomp-class}
-\title{Partially-observed Markov process}
+\keyword{internal}
+\title{Partially-observed Markov process class}
\description{
The class \code{pomp} encodes a partially-observed Markov process.
This page documents the structure of the class:
Modified: pkg/man/rmeasure-pomp.Rd
===================================================================
--- pkg/man/rmeasure-pomp.Rd 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/man/rmeasure-pomp.Rd 2010-04-29 19:54:45 UTC (rev 214)
@@ -3,6 +3,7 @@
\alias{rmeasure}
\alias{rmeasure,pomp-method}
\alias{rmeasure-pomp}
+\keyword{internal}
\title{Simulate the measurement model of a partially-observed Markov process}
\description{
The method \code{rmeasure} draws from the distribution of measurements given the state of the system.
Modified: pkg/man/rprocess-pomp.Rd
===================================================================
--- pkg/man/rprocess-pomp.Rd 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/man/rprocess-pomp.Rd 2010-04-29 19:54:45 UTC (rev 214)
@@ -3,7 +3,7 @@
\alias{rprocess}
\alias{rprocess,pomp-method}
\alias{rprocess-pomp}
-
+\keyword{internal}
\title{Simulate the process model of a partially-observed Markov process}
\description{
The method \code{rprocess} runs simulations of the the process-model portion of partially-observed Markov process.
Modified: pkg/man/skeleton-pomp.Rd
===================================================================
--- pkg/man/skeleton-pomp.Rd 2010-04-07 16:18:40 UTC (rev 213)
+++ pkg/man/skeleton-pomp.Rd 2010-04-29 19:54:45 UTC (rev 214)
@@ -3,7 +3,7 @@
\alias{skeleton}
\alias{skeleton,pomp-method}
\alias{skeleton-pomp}
-
+\keyword{internal}
\title{Evaluate the deterministic skeleton at the given points in state space.}
\description{
The method \code{skeleton} computes the deterministic skeleton.
More information about the pomp-commits
mailing list