[Pomp-commits] r215 - in pkg: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Apr 29 23:47:43 CEST 2010
Author: kingaa
Date: 2010-04-29 23:47:42 +0200 (Thu, 29 Apr 2010)
New Revision: 215
Modified:
pkg/R/compare.mif.R
pkg/R/mif.R
pkg/R/nlf-funcs.R
pkg/R/nlf-guts.R
pkg/R/nlf.R
pkg/R/pfilter-mif.R
pkg/R/pfilter.R
pkg/R/plot-pomp.R
pkg/R/pomp-methods.R
pkg/R/pomp.R
pkg/R/slice.R
pkg/R/trajectory-pomp.R
pkg/man/pfilter.Rd
Log:
- replace ':' construct with 'seq' and 'seq_len' throughout
- add 'seed' argument to 'pfilter'
Modified: pkg/R/compare.mif.R
===================================================================
--- pkg/R/compare.mif.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/compare.mif.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -27,7 +27,7 @@
time <- time(xx)
while (hi<nplots) {
hi <- min(low+n.per.page-1,nplots)
- for (i in low:hi) {
+ for (i in seq(from=low,to=hi,by=1)) {
n <- i-low+1
logplot <- if (plotnames[i]%in%lognames) "y" else ""
dat <- sapply(
@@ -75,7 +75,7 @@
iteration <- seq(0,xx at Nmif)
while (hi<nplots) {
hi <- min(low+n.per.page-1,nplots)
- for (i in low:hi) {
+ for (i in seq(from=low,to=hi,by=1)) {
n <- i-low+1
dat <- sapply(z,function(po,label) conv.rec(po,label),label=plotnames[i])
matplot(
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/mif.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -209,7 +209,7 @@
conv.rec=conv.rec
)
- for (n in seq(length=Nmif)) { # main loop
+ for (n in seq_len(Nmif)) { # main loop
## compute the cooled sigma
cool.sched <- try(
Modified: pkg/R/nlf-funcs.R
===================================================================
--- pkg/R/nlf-funcs.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/nlf-funcs.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -8,7 +8,7 @@
warning(" series length truncated to default in make.lags")
start <- max(lags)+1
temp <- matrix(0,ncol=xd*length(lags),nrow=n)
- for (k in 1:length(lags)) {
+ for (k in seq_len(length(lags))) {
a <- start-lags[k]
b <- a + n - 1
temp[,(1:xd)+(k-1)*xd] <- x[(a:b),]
@@ -48,7 +48,7 @@
X1 <- X-knots[1]
nknots <- length(knots)
if (nknots>1) {
- for (j in 2:nknots) {
+ for (j in seq(from=2,to=nknots,by=1)) {
X1 <- cbind(X1,X-knots[j])
}
}
@@ -68,7 +68,7 @@
Newey.West <- function(x,y,maxlag) {
w <- 1-(1:maxlag)/(maxlag+1)
out <- mean(x*y,na.rm=T)
- for(i in 1:maxlag) {
+ for (i in seq_len(maxlag)) {
out <- out+w[i]*mean(trimr(x,i,0)*trimr(y,0,i),na.rm=T)+w[i]*mean(trimr(y,i,0)*trimr(x,0,i),na.rm=T)
}
out
@@ -80,9 +80,9 @@
ncol.A <- ncol(A)
ncol.B <- ncol(B)
Tmat <- matrix(0,nrow(A),ncol.A*ncol.B)
- for(i in 1:ncol.A) {
+ for (i in seq_len(ncol.A)) {
start=(i-1)*ncol.B
- for(j in 1:ncol.B) {
+ for (j in seq_len(ncol.B)) {
Tmat[,start+j] <- A[,i]*B[,j]
}
}
Modified: pkg/R/nlf-guts.R
===================================================================
--- pkg/R/nlf-guts.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/nlf-guts.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -112,7 +112,7 @@
}
if (multivar) {
- for (jvar in 2:nvar) {
+ for (jvar in seq(from=2,to=nvar,by=1)) {
data.ts <- data.mat[jvar,]
model.ts <- model.mat[jvar,]
@@ -159,7 +159,7 @@
prediction.errors <- matrix(0,dim(data.pred)[1],nvar)
model.residuals <- matrix(0,dim(model.pred)[1],nvar)
- for (jvar in 1:nvar) {
+ for (jvar in seq_len(nvar)) {
model.lm <- lm(model.pred[,jvar]~rbfbasis.model-1)
model.residuals[,jvar] <- residuals(model.lm)
ck <- as.vector(coef(model.lm))
Modified: pkg/R/nlf.R
===================================================================
--- pkg/R/nlf.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/nlf.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -149,7 +149,7 @@
if (verbose) cat("h in NLF = ", h, "\n")
eps <- rep(h,nfitted)
- for (i in 1:nfitted) {
+ for (i in seq_len(nfitted)) {
Fvals <- rep(0,5)
Fvals[3] <- F0
guess <- fitted
@@ -182,7 +182,7 @@
if (verbose) cat("epsilon in NLF =",t(eps), "\n")
Imat <- matrix(0,npts,nfitted)
- for (i in 1:nfitted) {
+ for (i in seq_len(nfitted)) {
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,
@@ -210,8 +210,8 @@
Ihat[i,i] <- Newey.West(Imat[,i],Imat[,i],nlags)
}
- for (i in 1:(nfitted-1)) {
- for (j in (i+1):nfitted) {
+ for (i in seq_len(nfitted-1)) {
+ for (j in seq(from=i+1,to=nfitted,by=1)) {
guess.uu <- fitted
guess.uu[i] <- guess.uu[i]+eps[i]
guess.uu[j] <- guess.uu[j]+eps[j]
Modified: pkg/R/pfilter-mif.R
===================================================================
--- pkg/R/pfilter-mif.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/pfilter-mif.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -27,6 +27,7 @@
pred.var=pred.var,
filter.mean=filter.mean,
save.states=FALSE,
- ...)
+ ...
+ )
}
)
Modified: pkg/R/pfilter.R
===================================================================
--- pkg/R/pfilter.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/pfilter.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -6,8 +6,17 @@
pfilter.internal <- function (object, params, Np,
tol, max.fail,
pred.mean, pred.var, filter.mean,
- .rw.sd, verbose,
+ .rw.sd, seed, verbose,
save.states) {
+ if (missing(seed)) seed <- NULL
+ if (!is.null(seed)) {
+ if (!exists(".Random.seed",where=.GlobalEnv)) { # need to initialize the RNG
+ runif(1)
+ }
+ save.seed <- get(".Random.seed",pos=.GlobalEnv)
+ set.seed(seed)
+ }
+
if (missing(params)) {
params <- coef(object)
if (length(params)==0) {
@@ -98,7 +107,7 @@
)
else NULL
- for (nt in seq(length=ntimes)) {
+ for (nt in seq_len(ntimes)) {
## advance the state variables according to the process model
X <- try(
@@ -234,6 +243,12 @@
}
+ if (!is.null(seed)) {
+ assign(".Random.seed",save.seed,pos=.GlobalEnv)
+ seed <- save.seed
+ }
+
+
list(
pred.mean=pred.m,
pred.var=pred.v,
@@ -241,6 +256,7 @@
eff.sample.size=eff.sample.size,
cond.loglik=loglik,
states=xparticles,
+ seed=seed,
nfail=nfail,
loglik=sum(loglik)
)
@@ -256,6 +272,7 @@
pred.var = FALSE,
filter.mean = FALSE,
save.states = FALSE,
+ seed = NULL,
verbose = getOption("verbose"),
...) {
pfilter.internal(
@@ -268,6 +285,7 @@
pred.var=pred.var,
filter.mean=filter.mean,
save.states=save.states,
+ seed=seed,
verbose=verbose
)
}
Modified: pkg/R/plot-pomp.R
===================================================================
--- pkg/R/plot-pomp.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/plot-pomp.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -49,7 +49,7 @@
nr <- ceiling(nser/nc)
oldpar <- par(mar=mar,oma=oma,mfcol=c(nr,nc))
on.exit(par(oldpar))
- for (i in 1:nser) {
+ for (i in seq_len(nser)) {
plot.default(y=x[[i]],x=time,axes=FALSE,xlab="",ylab="",log=log,
col=col,bg=bg,pch=pch,ann=ann,type="n",...)
panel(y=x[[i]],x=time,col=col,bg=bg,pch=pch,type=type,...)
@@ -82,7 +82,7 @@
}
v1 <- 1
v2 <- min(v1+plots.per.page,length(vars))
- for (page in 1:n.page) {
+ for (page in seq_len(n.page)) {
vv <- vars[seq(from=v1,to=v2)]
plotpomp(
x=X[vv],
Modified: pkg/R/pomp-methods.R
===================================================================
--- pkg/R/pomp-methods.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/pomp-methods.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -83,7 +83,7 @@
)
if (ncol(ss)>length(tt)) tt <- c(tt0,tt)
nt <- c(object at t0,object at times)
- for (kt in seq(length=length(nt))) {
+ for (kt in seq_len(length(nt))) {
wr <- which(nt[kt]==tt)
if (length(wr)>0)
object at states[,kt] <- ss[,wr[1]]
Modified: pkg/R/pomp.R
===================================================================
--- pkg/R/pomp.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/pomp.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -290,7 +290,7 @@
nobs <- length(formulae)
if (nobs < 1)
stop("pomp error: to use ",sQuote("measurement.model")," you must provide at least one formula")
- for (k in 1:nobs) {
+ for (k in seq_len(nobs)) {
if (!inherits(formulae[[k]],"formula"))
stop("pomp error: ",sQuote("measurement.model")," takes formulae as arguments")
}
@@ -298,7 +298,7 @@
distrib <- lapply(formulae,function(x)as.character(x[[3]][[1]]))
ddistrib <- lapply(distrib,function(x)paste("d",x,sep=''))
rdistrib <- lapply(distrib,function(x)paste("r",x,sep=''))
- for (k in 1:nobs) {
+ for (k in seq_len(nobs)) {
res <- try(
match.fun(ddistrib[[k]]),
silent=TRUE
@@ -315,7 +315,7 @@
pred.args <- lapply(formulae,function(x)as.list(x[[3]][-1]))
dcalls <- vector(mode='list',length=nobs)
rcalls <- vector(mode='list',length=nobs)
- for (k in 1:nobs) {
+ for (k in seq_len(nobs)) {
dcalls[[k]] <- as.call(
c(
list(
@@ -341,7 +341,7 @@
list(
dmeasure = function (y, x, t, params, log, covars, ...) {
f <- 0
- for (k in 1:nobs) {
+ for (k in seq_len(nobs)) {
f <- f+eval(
dcalls[[k]],
envir=as.list(c(y,x,params,covars,t=t))
@@ -352,7 +352,7 @@
rmeasure = function (x, t, params, covars, ...) {
y <- numeric(length=nobs)
names(y) <- obsnames
- for (k in 1:nobs) {
+ for (k in seq_len(nobs)) {
y[k] <- eval(
rcalls[[k]],
envir=as.list(c(x,params,covars,t=t))
Modified: pkg/R/slice.R
===================================================================
--- pkg/R/slice.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/slice.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -10,7 +10,7 @@
ranges <- lapply(vars[varying],function(x)x[c(1,3)])
x <- as.data.frame(matrix(center.point,byrow=TRUE,ncol=nvars,nrow=n,dimnames=list(NULL,names(vars))))
y <- vector(mode="list",length=nranges)
- for (v in seq(length=nranges)) {
+ for (v in seq_len(nranges)) {
y[[v]] <- x
w <- varying[v]
y[[v]][[w]] <- seq(from=ranges[[v]][1],to=ranges[[v]][2],length=n)
@@ -25,8 +25,8 @@
y <- as.matrix(sobol(vars,n))
x <- as.matrix(do.call(expand.grid,prof))
z <- array(dim=c(nrow(x),nrow(y),ncol(x)+ncol(y)))
- for (j in 1:nrow(x)) {
- for (k in 1:nrow(y)) {
+ for (j in seq_len(nrow(x))) {
+ for (k in seq_len(nrow(y))) {
z[j,k,] <- c(x[j,],y[k,])
}
}
Modified: pkg/R/trajectory-pomp.R
===================================================================
--- pkg/R/trajectory-pomp.R 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/R/trajectory-pomp.R 2010-04-29 21:47:42 UTC (rev 215)
@@ -39,7 +39,7 @@
object at skeleton.type,
map={ # iterate the map
x[,,1] <- x0
- for (k in 2:length(times)) {
+ for (k in seq(from=2,to=length(times),by=1)) {
x[,,k] <- skeleton(
object,
x=x[,,k-1,drop=FALSE],
@@ -49,7 +49,7 @@
}
},
vectorfield={ # integrate the vectorfield
- for (j in 1:nrep) {
+ for (j in seq_len(nrep)) {
X <- try(
lsoda(
y=x0[,j],
Modified: pkg/man/pfilter.Rd
===================================================================
--- pkg/man/pfilter.Rd 2010-04-29 19:54:45 UTC (rev 214)
+++ pkg/man/pfilter.Rd 2010-04-29 21:47:42 UTC (rev 215)
@@ -7,13 +7,13 @@
\title{Particle filter}
\description{
Run a plain vanilla particle filter.
- Resampling is performed after each observation.
+ Resampling is performed at each observation.
}
\usage{
pfilter(object, \dots)
\S4method{pfilter}{pomp}(object, params, Np, tol = 1e-17,
max.fail = 0, pred.mean = FALSE, pred.var = FALSE,
- filter.mean = FALSE, save.states = FALSE,
+ filter.mean = FALSE, save.states = FALSE, seed = NULL,
verbose = getOption("verbose"), \dots)
\S4method{pfilter}{mif}(object, params, Np, tol = 1e-17,
max.fail = 0, pred.mean = FALSE, pred.var = FALSE,
@@ -25,7 +25,7 @@
}
\item{params}{
A \code{npars} x \code{Np} matrix containing the parameters corresponding to the initial state values in \code{xstart}.
- This must have a 'rownames' attribute.
+ This must have a \sQuote{rownames} attribute.
It is permissible to supply \code{params} as a named numeric vector, i.e., without a \code{dim} attribute.
In this case, all particles will inherit the same parameter values.
}
@@ -34,7 +34,7 @@
When \code{object} is of class \code{mif}, this is by default the same number of particles used in the \code{mif} iterations.
}
\item{tol}{
- positive numeric scalar; particles with log likelihood below \code{tol} are considered to be "lost".
+ positive numeric scalar; particles with log likelihood below \code{tol} are considered to be \dQuote{lost}.
A filtering failure occurs when, at some time point, all particles are lost.
When all particles are lost, the conditional log likelihood at that time point is set to be \code{log(tol)}.
}
@@ -54,6 +54,11 @@
\item{save.states}{
logical; if \code{TRUE}, the state-vector for each particle is saved and returned.
}
+ \item{seed}{
+ optional; an object specifying if and how the random number generator should be initialized (\sQuote{seeded}).
+ If \code{seed} is an integer, it is passed to \code{set.seed} prior to any simulation and is returned as the \dQuote{seed} element of the return list.
+ By default, the state of the random number generator is not changed and the value of \code{.Random.seed} on the call is stored in the \dQuote{seed} element of the return list.
+ }
\item{verbose}{
logical; if \code{TRUE}, progress information is reported as \code{pfilter} works.
}
@@ -84,6 +89,11 @@
An array with dimensions \code{nvars}-by-\code{Np}-by-\code{ntimes}.
In particular, \code{states[,i,t]} can be considered a sample from \eqn{f[X|y_{1:t}]}.
}
+ \item{seed}{
+ The state of the random number generator at the time \code{pfilter} was called.
+ If the argument \code{seed} was specified, this is a copy;
+ if not, this is the internal state of the random number generator at the time of call.
+ }
\item{nfail}{
The number of filtering failures encountered.
}
More information about the pomp-commits
mailing list