[Pomp-commits] r216 - in pkg: . R data inst inst/data-R man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 30 00:35:48 CEST 2010
Author: kingaa
Date: 2010-04-30 00:35:47 +0200 (Fri, 30 Apr 2010)
New Revision: 216
Modified:
pkg/DESCRIPTION
pkg/R/pfilter.R
pkg/R/trajmatch.R
pkg/data/euler.sir.rda
pkg/data/ou2.rda
pkg/data/rw2.rda
pkg/data/verhulst.rda
pkg/inst/ChangeLog
pkg/inst/data-R/ou2.R
pkg/man/trajmatch.Rd
pkg/src/skeleton.c
pkg/tests/examples.R
pkg/tests/logistic.Rout.save
pkg/tests/ou2-forecast.Rout.save
pkg/tests/ou2-kalman.Rout.save
pkg/tests/ou2-nlf.Rout.save
pkg/tests/ou2-simulate.Rout.save
pkg/tests/ou2-trajmatch.Rout.save
pkg/tests/rw2.Rout.save
pkg/tests/sir.Rout.save
Log:
- change the trajectory-matching algorithm so that it depends on 'trajectory': it is no longer to set the process noise to zero manually in order to use this method
- clean up error handling in 'pfilter'
- update .Rout.save files in 'tests'
- make 'skeleton' sensitive to variable names
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/DESCRIPTION 2010-04-29 22:35:47 UTC (rev 216)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.28-3
-Date: 2010-04-07
+Version: 0.28-4
+Date: 2010-04-29
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
Modified: pkg/R/pfilter.R
===================================================================
--- pkg/R/pfilter.R 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/R/pfilter.R 2010-04-29 22:35:47 UTC (rev 216)
@@ -20,7 +20,7 @@
if (missing(params)) {
params <- coef(object)
if (length(params)==0) {
- stop("pfilter error: ",sQuote("params")," must be supplied",call.=FALSE)
+ stop(sQuote("pfilter")," error: ",sQuote("params")," must be supplied",call.=FALSE)
}
}
if (missing(Np))
@@ -40,7 +40,7 @@
}
paramnames <- rownames(params)
if (is.null(paramnames))
- stop("pfilter error: ",sQuote("params")," must have rownames",call.=FALSE)
+ stop(sQuote("pfilter")," error: ",sQuote("params")," must have rownames",call.=FALSE)
x <- init.state(object,params=params)
statenames <- rownames(x)
@@ -59,9 +59,13 @@
if (random.walk) {
rw.names <- names(.rw.sd)
if (is.null(rw.names)||!is.numeric(.rw.sd))
- stop("pfilter error: ",sQuote(".rw.sd")," must be a named vector",call.=FALSE)
+ stop(sQuote("pfilter")," error: ",sQuote(".rw.sd")," must be a named vector",call.=FALSE)
if (any(!(rw.names%in%paramnames)))
- stop("pfilter error: the rownames of ",sQuote("params")," must include all of the names of ",sQuote(".rw.sd"),"",call.=FALSE)
+ stop(
+ sQuote("pfilter")," error: the rownames of ",
+ sQuote("params")," must include all of the names of ",
+ sQuote(".rw.sd"),"",call.=FALSE
+ )
sigma <- .rw.sd
} else {
rw.names <- character(0)
@@ -120,7 +124,7 @@
silent=FALSE
)
if (inherits(X,'try-error'))
- stop("pfilter error: process simulation error",call.=FALSE)
+ stop(sQuote("pfilter")," error: process simulation error",call.=FALSE)
x[,] <- X # ditch the third dimension
@@ -134,7 +138,7 @@
silent=FALSE
)
if (inherits(xx,'try-error')) {
- stop("pfilter error: error in prediction mean computation",call.=FALSE)
+ stop(sQuote("pfilter")," error: error in prediction mean computation",call.=FALSE)
} else {
pred.m[,nt] <- xx
}
@@ -145,7 +149,7 @@
problem.indices <- unique(which(!is.finite(x),arr.ind=TRUE)[,1])
if (length(problem.indices)>0) {
stop(
- "pfilter error: non-finite state variable(s): ",
+ sQuote("pfilter")," error: non-finite state variable(s): ",
paste(rownames(x)[problem.indices],collapse=', '),
call.=FALSE
)
@@ -154,7 +158,7 @@
problem.indices <- unique(which(!is.finite(params[rw.names,,drop=FALSE]),arr.ind=TRUE)[,1])
if (length(problem.indices)>0) {
stop(
- "pfilter error: non-finite parameter(s): ",
+ sQuote("pfilter")," error: non-finite parameter(s): ",
paste(rw.names[problem.indices],collapse=', '),
call.=FALSE
)
@@ -168,7 +172,7 @@
silent=FALSE
)
if (inherits(xx,'try-error')) {
- stop("pfilter error: error in prediction variance computation",call.=FALSE)
+ stop(sQuote("pfilter")," error: error in prediction variance computation",call.=FALSE)
} else {
pred.v[,nt] <- xx
}
@@ -186,10 +190,10 @@
silent=FALSE
)
if (inherits(weights,'try-error'))
- stop("pfilter error: error in calculation of weights",call.=FALSE)
+ stop(sQuote("pfilter")," error: error in calculation of weights",call.=FALSE)
if (any(is.na(weights))) {
## problem.indices <- which(is.na(weights))
- stop("pfilter error: dmeasure returns NA",call.=FALSE)
+ stop(sQuote("pfilter")," error: ",sQuote("dmeasure")," returns NA",call.=FALSE)
}
## test for failure to filter
@@ -201,7 +205,7 @@
message("filtering failure at time t = ",times[nt+1])
nfail <- nfail+1
if (nfail > max.fail)
- stop('pfilter error: too many filtering failures',call.=FALSE)
+ stop(sQuote("pfilter")," error: too many filtering failures",call.=FALSE)
loglik[nt] <- log(tol) # worst log-likelihood
weights <- rep(1/Np,Np)
eff.sample.size[nt] <- 0
Modified: pkg/R/trajmatch.R
===================================================================
--- pkg/R/trajmatch.R 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/R/trajmatch.R 2010-04-29 22:35:47 UTC (rev 216)
@@ -17,7 +17,8 @@
fn=function (x) {
p <- start
p[par.est] <- x
- x <- simulate(object,nsim=1,params=p,states=TRUE)[,,-1,drop=FALSE]
+## x <- simulate(object,nsim=1,params=p,states=TRUE)[,,-1,drop=FALSE]
+ x <- trajectory(object,params=p)[,,-1,drop=FALSE]
d <- dmeasure(
object,
y=data.array(object),
Modified: pkg/data/euler.sir.rda
===================================================================
(Binary files differ)
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/rw2.rda
===================================================================
(Binary files differ)
Modified: pkg/data/verhulst.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/inst/ChangeLog 2010-04-29 22:35:47 UTC (rev 216)
@@ -1,5 +1,26 @@
+2010-04-29 kingaa
+
+ * [r215] R/compare.mif.R, R/mif.R, R/nlf-funcs.R, R/nlf-guts.R,
+ R/nlf.R, R/pfilter-mif.R, R/pfilter.R, R/plot-pomp.R,
+ R/pomp-methods.R, R/pomp.R, R/slice.R, R/trajectory-pomp.R,
+ man/pfilter.Rd: - replace ':' construct with 'seq' and 'seq_len'
+ throughout
+ - add 'seed' argument to 'pfilter'
+ * [r214] R/nlf-funcs.R, R/nlf-lql.R, R/nlf-objfun.R, R/nlf.R,
+ man/dmeasure-pomp.Rd, man/dprocess-pomp.Rd,
+ man/init.state-pomp.Rd, man/mif-class.Rd, man/nlf.Rd,
+ man/particles-mif.Rd, man/pomp-class.Rd, man/rmeasure-pomp.Rd,
+ man/rprocess-pomp.Rd, man/skeleton-pomp.Rd: - 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
+
2010-04-07 kingaa
+ * [r213] DESCRIPTION, inst/ChangeLog,
+ inst/doc/advanced_topics_in_pomp.pdf, inst/doc/intro_to_pomp.pdf:
+ - version 0.28-3
* [r212] man/LondonYorke.Rd: - improve documentation slightly
2010-04-06 kingaa
Modified: pkg/inst/data-R/ou2.R
===================================================================
--- pkg/inst/data-R/ou2.R 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/inst/data-R/ou2.R 2010-04-29 22:35:47 UTC (rev 216)
@@ -56,6 +56,17 @@
},
dmeasure = "normal_dmeasure",
rmeasure = "normal_rmeasure",
+ skeleton.map = function (x, t, params, ...) {
+ with(
+ as.list(c(x,params)),
+ {
+ c(
+ x1=alpha.1*x1+alpha.3*x2,
+ x2=alpha.2*x1+alpha.4*x2
+ )
+ }
+ )
+ },
paramnames = c(
"alpha.1","alpha.2","alpha.3","alpha.4",
"sigma.1","sigma.2","sigma.3",
Modified: pkg/man/trajmatch.Rd
===================================================================
--- pkg/man/trajmatch.Rd 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/man/trajmatch.Rd 2010-04-29 22:35:47 UTC (rev 216)
@@ -23,29 +23,31 @@
}
\details{
Trajectory matching is accomplished using \code{\link{optim}}.
- It is assumed that the process model is deterministic.
+ The \code{\link{trajectory}} method is used for this, which in turn uses the \code{skeleton} slot of the \code{pomp} object.
}
\examples{
data(ou2)
true.p <- c(
- alpha.1=0.9,alpha.2=0,alpha.3=0,alpha.4=0.99,
- sigma.1=1,sigma.2=0,sigma.3=2,
+ alpha.1=0.9,alpha.2=0,alpha.3=-0.4,alpha.4=0.99,
+ sigma.1=2,sigma.2=0.1,sigma.3=2,
tau=1,
x1.0=50,x2.0=-50
)
simdata <- simulate(ou2,nsim=1,params=true.p,seed=43553)
guess.p <- true.p
- guess.p[grep('sigma',names(guess.p))] <- 0 ## make the process model deterministic
res <- traj.match(
simdata,
start=guess.p,
- est=c('alpha.1','alpha.4','x1.0','x2.0','tau'),
+ est=c('alpha.1','alpha.3','alpha.4','x1.0','x2.0','tau'),
maxit=2000,
+ method="Nelder-Mead",
reltol=1e-8
)
+
}
\seealso{
\code{\link{optim}},
+ \code{\link{trajectory}},
\code{\link{pomp}},
\code{\link{pomp-class}}
}
Modified: pkg/src/skeleton.c
===================================================================
--- pkg/src/skeleton.c 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/src/skeleton.c 2010-04-29 22:35:47 UTC (rev 216)
@@ -58,6 +58,8 @@
static int _pomp_skel_npar; // number of parameters
static SEXP _pomp_skel_envir; // function's environment
static SEXP _pomp_skel_fcall; // function call
+static SEXP _pomp_skel_vnames; // names of state variables
+//static int *_pomp_skel_vindex; // indices of state variables
#define XVEC (_pomp_skel_Xvec)
#define PVEC (_pomp_skel_Pvec)
@@ -67,6 +69,8 @@
#define NPAR (_pomp_skel_npar)
#define RHO (_pomp_skel_envir)
#define FCALL (_pomp_skel_fcall)
+#define VNAMES (_pomp_skel_vnames)
+//#define VINDEX (_pomp_skel_vindex)
// this is the vectorfield that is evaluated when the user supplies an R function
// (and not a native routine)
@@ -78,7 +82,9 @@
int nprotect = 0;
int k;
double *xp;
- SEXP ans;
+ SEXP ans, nm, idx;
+ int use_names, *op;
+
xp = REAL(XVEC);
for (k = 0; k < NVAR; k++) xp[k] = x[k];
xp = REAL(PVEC);
@@ -88,12 +94,29 @@
xp = REAL(TIME);
xp[0] = t;
PROTECT(ans = eval(FCALL,RHO)); nprotect++; // evaluate the call
+
if (LENGTH(ans)!=NVAR) {
UNPROTECT(nprotect);
error("skeleton error: user 'skeleton' must return a vector of length %d",NVAR);
}
+
+ PROTECT(nm = GET_NAMES(ans)); nprotect++;
+ use_names = !isNull(nm);
+ if (use_names) { // match names against names from data slot
+ PROTECT(idx = matchnames(VNAMES,nm)); nprotect++;
+ op = INTEGER(idx);
+ // for (k = 0; k < NVAR; k++) VINDEX[k] = op[k];
+ // } else {
+ // VINDEX = 0;
+ }
+
xp = REAL(AS_NUMERIC(ans));
- for (k = 0; k < NVAR; k++) f[k] = xp[k];
+ if (use_names) {
+ for (k = 0; k < NVAR; k++) f[op[k]] = xp[k];
+ } else {
+ for (k = 0; k < NVAR; k++) f[k] = xp[k];
+ }
+
UNPROTECT(nprotect);
}
@@ -109,7 +132,7 @@
SEXP tcovar, covar;
SEXP statenames, paramnames, covarnames;
SEXP sindex, pindex, cindex;
- SEXP Xnames, Pnames, Cnames;
+ SEXP Pnames, Cnames;
pomp_vectorfield_map *ff = NULL;
int k;
@@ -150,7 +173,7 @@
ncovars = LENGTH(covarnames);
dim = INTEGER(GET_DIM(covar)); covlen = dim[0]; covdim = dim[1];
- PROTECT(Xnames = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
+ PROTECT(VNAMES = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
PROTECT(Cnames = GET_COLNAMES(GET_DIMNAMES(covar))); nprotect++;
@@ -174,7 +197,7 @@
for (k = 0; k < nvars; k++) xp[k] = 0.0;
PROTECT(PVEC = NEW_NUMERIC(npars)); nprotect++; // for internal use
PROTECT(CVEC = NEW_NUMERIC(covdim)); nprotect++; // for internal use
- SET_NAMES(XVEC,Xnames); // make sure the names attribute is copied
+ SET_NAMES(XVEC,VNAMES); // make sure the names attribute is copied
SET_NAMES(PVEC,Pnames); // make sure the names attribute is copied
SET_NAMES(CVEC,Cnames); // make sure the names attribute is copied
// set up the function call
@@ -192,7 +215,7 @@
fdim[0] = nvars; fdim[1] = nreps; fdim[2] = ntimes;
PROTECT(F = makearray(3,fdim)); nprotect++;
- setrownames(F,Xnames,3);
+ setrownames(F,VNAMES,3);
xp = REAL(F);
for (k = 0; k < nvars*nreps*ntimes; k++) xp[k] = 0.0;
Modified: pkg/tests/examples.R
===================================================================
--- pkg/tests/examples.R 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/tests/examples.R 2010-04-29 22:35:47 UTC (rev 216)
@@ -5,5 +5,5 @@
examples <- list.files(path=system.file("examples",package="pomp"),pattern=".\\.R$",full.names=TRUE)
for (e in examples) {
- source(e,local=TRUE)
+ source(e,local=TRUE,echo=TRUE)
}
Modified: pkg/tests/logistic.Rout.save
===================================================================
--- pkg/tests/logistic.Rout.save 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/tests/logistic.Rout.save 2010-04-29 22:35:47 UTC (rev 216)
@@ -1,6 +1,6 @@
-R version 2.8.1 (2008-12-22)
-Copyright (C) 2008 The R Foundation for Statistical Computing
+R version 2.10.1 (2009-12-14)
+Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
Modified: pkg/tests/ou2-forecast.Rout.save
===================================================================
--- pkg/tests/ou2-forecast.Rout.save 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/tests/ou2-forecast.Rout.save 2010-04-29 22:35:47 UTC (rev 216)
@@ -1,6 +1,6 @@
-R version 2.11.0 Under development (unstable) (2010-02-25 r51179)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.10.1 (2009-12-14)
+Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
Modified: pkg/tests/ou2-kalman.Rout.save
===================================================================
--- pkg/tests/ou2-kalman.Rout.save 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/tests/ou2-kalman.Rout.save 2010-04-29 22:35:47 UTC (rev 216)
@@ -1,6 +1,6 @@
-R version 2.8.1 (2008-12-22)
-Copyright (C) 2008 The R Foundation for Statistical Computing
+R version 2.10.1 (2009-12-14)
+Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -104,7 +104,7 @@
> kalm.fit1 <- optim(p.guess,kalman,object=ou2,params=p.truth,hessian=T)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 37.90149 secs
+Time difference of 17.08114 secs
> tic <- Sys.time()
> print(-kalm.fit1$value,digits=4)
[1] -411.1
Modified: pkg/tests/ou2-nlf.Rout.save
===================================================================
--- pkg/tests/ou2-nlf.Rout.save 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/tests/ou2-nlf.Rout.save 2010-04-29 22:35:47 UTC (rev 216)
@@ -1,5 +1,5 @@
-R version 2.9.1 (2009-06-26)
+R version 2.10.1 (2009-12-14)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
@@ -42,7 +42,7 @@
+ )
>
> print(m1)
-[1] -390.961
+[1] -391.5853
>
> m1 <- nlf(
+ object=po,
@@ -58,56 +58,62 @@
+ lql.frac = 0.025
+ )
Nelder-Mead direct search function minimizer
-function value for initial parameters = 390.960996
- Scaled convergence tolerance is 5.82577e-06
+function value for initial parameters = 391.585287
+ Scaled convergence tolerance is 5.83508e-06
Stepsize computed as 0.020000
-BUILD 3 391.170142 390.960996
-EXTENSION 5 391.121379 390.848054
-EXTENSION 7 390.960996 390.495673
-LO-REDUCTION 9 390.848054 390.495673
-EXTENSION 11 390.509255 390.032148
-EXTENSION 13 390.495673 389.677144
-REFLECTION 15 390.032148 389.601518
-EXTENSION 17 389.677144 389.069152
-LO-REDUCTION 19 389.601518 389.069152
-LO-REDUCTION 21 389.316744 389.069152
-REFLECTION 23 389.097526 388.964593
-LO-REDUCTION 25 389.069152 388.964593
-LO-REDUCTION 27 388.977691 388.959018
-HI-REDUCTION 29 388.964593 388.955938
-LO-REDUCTION 31 388.959018 388.954600
-HI-REDUCTION 33 388.955938 388.954600
-HI-REDUCTION 35 388.955009 388.953977
-LO-REDUCTION 37 388.954600 388.953960
-HI-REDUCTION 39 388.954009 388.953960
-HI-REDUCTION 41 388.953977 388.953875
-HI-REDUCTION 43 388.953960 388.953849
-LO-REDUCTION 45 388.953875 388.953849
-HI-REDUCTION 47 388.953868 388.953849
+BUILD 3 391.822251 391.585287
+EXTENSION 5 391.776151 391.437121
+EXTENSION 7 391.585287 390.985005
+EXTENSION 9 391.437121 390.705228
+EXTENSION 11 390.985005 389.914394
+LO-REDUCTION 13 390.705228 389.914394
+REFLECTION 15 390.162389 389.691267
+EXTENSION 17 389.914394 389.155652
+LO-REDUCTION 19 389.691267 389.155652
+EXTENSION 21 389.331663 388.843884
+REFLECTION 23 389.155652 388.683945
+LO-REDUCTION 25 388.843884 388.683945
+LO-REDUCTION 27 388.824168 388.683945
+REFLECTION 29 388.701174 388.673612
+HI-REDUCTION 31 388.683945 388.665874
+REFLECTION 33 388.673612 388.656432
+LO-REDUCTION 35 388.665874 388.656054
+LO-REDUCTION 37 388.656432 388.654115
+HI-REDUCTION 39 388.656054 388.651779
+HI-REDUCTION 41 388.654115 388.650952
+LO-REDUCTION 43 388.651779 388.650931
+HI-REDUCTION 45 388.650952 388.650792
+HI-REDUCTION 47 388.650931 388.650773
+HI-REDUCTION 49 388.650792 388.650754
+HI-REDUCTION 51 388.650773 388.650704
+LO-REDUCTION 53 388.650754 388.650704
+HI-REDUCTION 55 388.650718 388.650704
+HI-REDUCTION 57 388.650715 388.650704
+REFLECTION 59 388.650710 388.650703
Exiting from Nelder Mead minimizer
- 49 function evaluations used
+ 61 function evaluations used
h in NLF = 0.1
-epsilon in NLF = 0.3120962 0.3004542
-Fitted param 1 -3.984644 -3.984644 up in 'nlf'
-Fitted param 1 -3.985982 down in 'NLF'
-Fitted param 2 -3.993283 -3.993283 up in 'nlf'
-Fitted param 2 -3.993728 down in 'NLF'
+epsilon in NLF = 0.3450102 0.291259
+Fitted param 1 -3.986327 -3.986327 up in 'nlf'
+Fitted param 1 -3.999289 down in 'NLF'
+Fitted param 2 -3.990065 -3.990065 up in 'nlf'
+Fitted param 2 -3.991384 down in 'NLF'
>
> se <- p.truth
> se[] <- NA
> se[names(m1$se)] <- m1$se
> print(cbind(truth=p.truth,fit=m1$params,se=se),digits=3)
- truth fit se
-alpha.1 0.1 -0.1539 0.114
-alpha.2 0.0 0.0000 NA
-alpha.3 0.0 0.0000 NA
-alpha.4 0.2 0.0159 0.110
-sigma.1 1.0 1.0000 NA
-sigma.2 0.0 0.0000 NA
-sigma.3 2.0 2.0000 NA
-tau 1.0 1.0000 NA
-x1.0 0.0 0.0000 NA
-x2.0 0.0 0.0000 NA
+ truth fit se
+alpha.1 0.1 -0.2634 0.0894
+alpha.2 0.0 0.0000 NA
+alpha.3 0.0 0.0000 NA
+alpha.4 0.2 0.0389 0.1060
+sigma.1 1.0 1.0000 NA
+sigma.2 0.0 0.0000 NA
+sigma.3 2.0 2.0000 NA
+tau 1.0 1.0000 NA
+x1.0 0.0 0.0000 NA
+x2.0 0.0 0.0000 NA
>
> po <- simulate(po,times=(1:1000))
> m2 <- nlf(
@@ -124,45 +130,47 @@
+ lql.frac = 0.025
+ )
Nelder-Mead direct search function minimizer
-function value for initial parameters = 3931.875269
- Scaled convergence tolerance is 5.85895e-05
+function value for initial parameters = 4016.057678
+ Scaled convergence tolerance is 5.98439e-05
Stepsize computed as 0.020000
-BUILD 3 3931.875269 3931.417243
-EXTENSION 5 3931.541801 3930.827423
-EXTENSION 7 3931.417243 3930.572713
-REFLECTION 9 3930.827423 3930.228048
-LO-REDUCTION 11 3930.572713 3930.228048
-LO-REDUCTION 13 3930.344970 3930.215357
-HI-REDUCTION 15 3930.238133 3930.215357
-HI-REDUCTION 17 3930.228048 3930.213344
-HI-REDUCTION 19 3930.215357 3930.209340
-HI-REDUCTION 21 3930.213344 3930.208245
-LO-REDUCTION 23 3930.209340 3930.207339
-HI-REDUCTION 25 3930.208245 3930.207190
-HI-REDUCTION 27 3930.207339 3930.206898
-HI-REDUCTION 29 3930.207190 3930.206898
+BUILD 3 4016.416157 4016.057678
+REFLECTION 5 4016.318308 4016.057663
+EXTENSION 7 4016.057678 4015.676533
+LO-REDUCTION 9 4016.057663 4015.676533
+REFLECTION 11 4015.681034 4015.580294
+HI-REDUCTION 13 4015.676533 4015.580294
+REFLECTION 15 4015.581259 4015.556202
+HI-REDUCTION 17 4015.580294 4015.537804
+LO-REDUCTION 19 4015.556202 4015.537804
+HI-REDUCTION 21 4015.539048 4015.537770
+HI-REDUCTION 23 4015.537804 4015.533394
+HI-REDUCTION 25 4015.537770 4015.532439
+LO-REDUCTION 27 4015.533394 4015.532436
+HI-REDUCTION 29 4015.532439 4015.532223
+HI-REDUCTION 31 4015.532436 4015.532160
+HI-REDUCTION 33 4015.532223 4015.532159
Exiting from Nelder Mead minimizer
- 31 function evaluations used
+ 35 function evaluations used
h in NLF = 0.1
-epsilon in NLF = 0.4209784 0.2679104
-Fitted param 1 -3.976024 -3.976024 up in 'nlf'
-Fitted param 1 -3.963879 down in 'NLF'
-Fitted param 2 -3.967971 -3.967971 up in 'nlf'
-Fitted param 2 -3.966288 down in 'NLF'
+epsilon in NLF = 0.4611727 0.2676601
+Fitted param 1 -4.057614 -4.057614 up in 'nlf'
+Fitted param 1 -4.054204 down in 'NLF'
+Fitted param 2 -4.052918 -4.052918 up in 'nlf'
+Fitted param 2 -4.051148 down in 'NLF'
>
> se <- p.truth
> se[] <- NA
> se[names(m2$se)] <- m2$se
> print(cbind(truth=p.truth,fit=m2$params,se=se),digits=3)
- truth fit se
-alpha.1 0.1 0.196 0.0541
-alpha.2 0.0 0.000 NA
-alpha.3 0.0 0.000 NA
-alpha.4 0.2 0.234 0.0372
-sigma.1 1.0 1.000 NA
-sigma.2 0.0 0.000 NA
-sigma.3 2.0 2.000 NA
-tau 1.0 1.000 NA
-x1.0 0.0 0.000 NA
-x2.0 0.0 0.000 NA
+ truth fit se
+alpha.1 0.1 0.0349 0.0662
+alpha.2 0.0 0.0000 NA
+alpha.3 0.0 0.0000 NA
+alpha.4 0.2 0.1916 0.0523
+sigma.1 1.0 1.0000 NA
+sigma.2 0.0 0.0000 NA
+sigma.3 2.0 2.0000 NA
+tau 1.0 1.0000 NA
+x1.0 0.0 0.0000 NA
+x2.0 0.0 0.0000 NA
>
Modified: pkg/tests/ou2-simulate.Rout.save
===================================================================
--- pkg/tests/ou2-simulate.Rout.save 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/tests/ou2-simulate.Rout.save 2010-04-29 22:35:47 UTC (rev 216)
@@ -1,6 +1,6 @@
-R version 2.8.1 (2008-12-22)
-Copyright (C) 2008 The R Foundation for Statistical Computing
+R version 2.10.1 (2009-12-14)
+Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -34,7 +34,7 @@
> ou2.sim <- simulate(ou2,params=p,nsim=100,seed=32043858)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.3518229 secs
+Time difference of 0.1700060 secs
>
> coef(ou2,c('x1.0','x2.0')) <- c(-50,50)
>
Modified: pkg/tests/ou2-trajmatch.Rout.save
===================================================================
--- pkg/tests/ou2-trajmatch.Rout.save 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/tests/ou2-trajmatch.Rout.save 2010-04-29 22:35:47 UTC (rev 216)
@@ -1,5 +1,5 @@
-R version 2.9.1 (2009-06-26)
+R version 2.10.1 (2009-12-14)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Modified: pkg/tests/rw2.Rout.save
===================================================================
--- pkg/tests/rw2.Rout.save 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/tests/rw2.Rout.save 2010-04-29 22:35:47 UTC (rev 216)
@@ -1,5 +1,5 @@
-R version 2.9.1 (2009-06-26)
+R version 2.10.1 (2009-12-14)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
@@ -245,18 +245,18 @@
{
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)))
}
y
}
-<environment: 0x27c8440>
+<environment: 0x189dbb0>
measurement model density, 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)))
}
@@ -264,7 +264,7 @@
f
else exp(f)
}
-<environment: 0x27c8440>
+<environment: 0x189dbb0>
initializer =
function (params, t0, ...)
{
Modified: pkg/tests/sir.Rout.save
===================================================================
--- pkg/tests/sir.Rout.save 2010-04-29 21:47:42 UTC (rev 215)
+++ pkg/tests/sir.Rout.save 2010-04-29 22:35:47 UTC (rev 216)
@@ -1,6 +1,6 @@
-R version 2.11.0 Under development (unstable) (2010-02-25 r51179)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.10.1 (2009-12-14)
+Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -417,18 +417,18 @@
{
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)))
}
y
}
-<environment: 0x2c7b448>
+<environment: 0x325f890>
measurement model density, 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)))
}
@@ -436,7 +436,7 @@
f
else exp(f)
}
-<environment: 0x2c7b448>
+<environment: 0x325f890>
initializer =
function (params, t0, ...)
{
@@ -507,7 +507,7 @@
> x <- simulate(po,params=log(params),nsim=3)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.009282 secs
+Time difference of 2.060826 secs
>
> pdf(file='sir.pdf')
>
@@ -524,7 +524,7 @@
> X3 <- trajectory(po,params=log(params),times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 4.513309 secs
+Time difference of 4.699151 secs
> plot(t3,X3['I',1,],type='l')
>
> f1 <- dprocess(
@@ -584,7 +584,7 @@
> x <- simulate(po,nsim=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 2.146718 secs
+Time difference of 2.14945 secs
> plot(x[[1]],variables=c("S","I","R","cases","W"))
>
> t3 <- seq(0,20,by=1/52)
@@ -592,7 +592,7 @@
> X4 <- trajectory(po,times=t3,hmax=1/52)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 0.7040031 secs
+Time difference of 0.734072 secs
> plot(t3,X4['I',1,],type='l')
>
> g2 <- dmeasure(
More information about the pomp-commits
mailing list