[Pomp-commits] r558 - in pkg: . R inst man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Sep 5 17:37:23 CEST 2011
Author: kingaa
Date: 2011-09-05 17:37:23 +0200 (Mon, 05 Sep 2011)
New Revision: 558
Added:
pkg/tests/ou2-icfit.R
pkg/tests/ou2-icfit.Rout.save
Removed:
pkg/tests/sir-icfit.R
pkg/tests/sir-icfit.Rout.save
Modified:
pkg/DESCRIPTION
pkg/R/mif.R
pkg/inst/NEWS
pkg/man/mif.Rd
Log:
- 'mif' is modified to allow estimation of IVPs only
- update the 'mif' help page to document this
- add codes to test this in 'tests/ou2-icfit'
- remove the 'sir-icfit' test codes for the moment
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2011-09-03 13:54:05 UTC (rev 557)
+++ pkg/DESCRIPTION 2011-09-05 15:37:23 UTC (rev 558)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.39-3
-Date: 2011-09-03
+Date: 2011-09-05
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>
URL: http://pomp.r-forge.r-project.org
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2011-09-03 13:54:05 UTC (rev 557)
+++ pkg/R/mif.R 2011-09-05 15:37:23 UTC (rev 558)
@@ -65,9 +65,6 @@
if (missing(pars))
stop("mif error: ",sQuote("pars")," must be specified",call.=FALSE)
- if (length(pars)==0)
- stop("mif error: at least one ordinary (non-IVP) parameter must be estimated",call.=FALSE)
-
if (missing(ivps))
stop("mif error: ",sQuote("ivps")," must be specified",call.=FALSE)
@@ -134,6 +131,23 @@
ic.lag <- as.integer(ic.lag)
if ((length(ic.lag)!=1)||(ic.lag < 1))
stop("mif error: ",sQuote("ic.lag")," must be a positive integer",call.=FALSE)
+ if (ic.lag>ntimes) {
+ warning(
+ "mif warning: ",sQuote("ic.lag")," = ",ic.lag," > ",ntimes,
+ " = length(time(",sQuote("object"),"))",
+ " is nonsensical. Setting ",sQuote("ic.lag")," = ",ntimes,".",
+ call.=FALSE
+ )
+ ic.lag <- length(time(object))
+ }
+ if ((length(pars)==0)&&(ic.lag<length(time(object)))) {
+ warning(
+ "mif warning: only IVPs are to be estimated, yet ",sQuote("ic.lag")," = ",ic.lag,
+ " < ",ntimes," = length(time(",sQuote("object"),")),",
+ " so unnecessary work is to be done.",
+ call.=FALSE
+ )
+ }
if (missing(cooling.factor))
stop("mif error: ",sQuote("cooling.factor")," must be specified",call.=FALSE)
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2011-09-03 13:54:05 UTC (rev 557)
+++ pkg/inst/NEWS 2011-09-05 15:37:23 UTC (rev 558)
@@ -5,6 +5,8 @@
o Inside 'init.state', a check has been added to make sure that the user's initializer function returns a vector of uniform size, as per the algorithms' assumptions.
Thanks to Micaela Martinez-Bakker for catching this.
+ o 'mif' has been modified so as to allow estimation of initial-value parameters alone via fixed-lag smoothing.
+
0.39-2
o When 'po' is a 'pomp'-class object with covariates, 'as.data.frame(po)' and 'as(po,"data.frame")' now contain columns for the covariates.
The latter are interpolated, if necessary, at the observation times.
Modified: pkg/man/mif.Rd
===================================================================
--- pkg/man/mif.Rd 2011-09-03 13:54:05 UTC (rev 557)
+++ pkg/man/mif.Rd 2011-09-05 15:37:23 UTC (rev 558)
@@ -53,8 +53,9 @@
Leaving \code{pars} unspecified is equivalent to setting it equal to the names of all parameters with a positive value of \code{rw.sd} that are not \code{ivps}.
}
\item{ivps}{
- optional character vector naming the initial-value parameters to be estimated.
+ optional character vector naming the initial-value parameters (IVPs) to be estimated.
Every parameter named in \code{ivps} must have a positive random-walk standard deviation specified in \code{rw.sd}.
+ If \code{pars} is empty, i.e., only IVPs are to be estimated, see below \dQuote{"Using MIF to estimate initial-value parameters only"}.
}
\item{particles}{
Function of prototype \code{particles(Np,center,sd,...)} which sets up the starting particle matrix by drawing a sample of size \code{Np} from the starting particle distribution centered at \code{center} and of width \code{sd}.
@@ -74,7 +75,7 @@
\item{Np}{
the number of particles to use in filtering.
This may be specified as a single positive integer, in which case the same number of particles will be used at each timestep.
- Alternatively, if one wishes the number of particles to vary across timesteps, one may specify \code{Np} either as a vector of positive integers (of length \code{length(time(object,t0=TRUE))}) or as a function taking a positive integer argument.
+ Alternatively, if one wishes the number of particles to vary across timestep, one may specify \code{Np} either as a vector of positive integers (of length \code{length(time(object,t0=TRUE))}) or as a function taking a positive integer argument.
In the latter case, \code{Np(k)} must be a single positive integer, representing the number of particles to be used at the \code{k}-th timestep:
\code{Np(0)} is the number of particles to use going from \code{timezero(object)} to \code{time(object)[1]},
\code{Np(1)}, from \code{timezero(object)} to \code{time(object)[1]},
@@ -85,6 +86,8 @@
a positive integer;
the timepoint for fixed-lag smoothing of initial-value parameters.
The \code{mif} update for initial-value parameters consists of replacing them by their filtering mean at time \code{times[ic.lag]}, where \code{times=time(object)}.
+ It makes no sense to set \code{ic.lag>length(times)};
+ if it is so set, \code{ic.lag} is set to \code{length(times)} with a warning.
}
\item{var.factor}{
a positive number;
@@ -130,11 +133,18 @@
If one does specify additional arguments, these will override the defaults.
}
\section{Continuing MIF Iterations}{
- One can continue a series of MIF iterations from where one left off using the \code{continue} method.
+ One can resume a series of MIF iterations from where one left off using the \code{continue} method.
A call to \code{mif} to perform \code{Nmif=m} iterations followed by a call to \code{continue} to perform \code{Nmif=n} iterations will produce precisely the same effect as a single call to \code{mif} to perform \code{Nmif=m+n} iterations.
By default, all the algorithmic parameters are the same as used in the original call to \code{mif}.
Additional arguments will override the defaults.
}
+\section{Using MIF to estimate initial-value parameters only}{
+ One can use MIF's fixed-lag smoothing to estimate only initial value parameters (IVPs).
+ In this case, \code{pars} is left empty and the IVPs to be estimated are named in \code{ivps}.
+ If \code{theta} is the current parameter vector, then at each MIF iteration, \code{Np} particles are drawn from a distribution centered at \code{theta} and with width proportional to \code{var.factor*rw.sd}, a particle filtering operation is performed, and \code{theta} is replaced by the filtering mean at \code{time(object)[ic.lag]}.
+ Note the implication that, when \code{mif} is used in this way on a time series any longer than \code{ic.lag}, unnecessary work is done.
+ If the time series in \code{object} is longer than \code{ic.lag}, consider replacing \code{object} with \code{window(object,end=ic.lag)}.
+}
% \section{Likelihood profiles using MIF}{
% The function \code{mif.profile.design} sets up a MIF likelihood-profile calculation.
% It returns a list of lists, each one of which contains a \code{mif} object.
Added: pkg/tests/ou2-icfit.R
===================================================================
--- pkg/tests/ou2-icfit.R (rev 0)
+++ pkg/tests/ou2-icfit.R 2011-09-05 15:37:23 UTC (rev 558)
@@ -0,0 +1,64 @@
+library(pomp)
+
+set.seed(921625222L)
+
+data(ou2)
+
+ics <- c("x1.0","x2.0") # names of the initial condition parameters
+
+p <- coef(ou2)
+p[ics] <- c(4,-4)
+
+fit <- mif(
+ ou2,
+ start=p,
+ Nmif=2,
+ ic.lag=1000,
+ var.factor=4,
+ ivps=ics,
+ rw.sd=c(
+ x1.0=1,x2.0=1
+ ),
+ Np=1000,
+ cooling.factor=1,
+ max.fail=10
+ )
+
+fit <- mif(
+ ou2,
+ start=p,
+ Nmif=2,
+ ic.lag=10,
+ var.factor=4,
+ ivps=ics,
+ rw.sd=c(
+ x1.0=1,x2.0=1
+ ),
+ Np=1000,
+ cooling.factor=1,
+ max.fail=10
+ )
+
+fit <- mif(
+ window(ou2,end=10),
+ start=p,
+ Nmif=500,
+ ic.lag=10,
+ var.factor=4,
+ ivps=ics,
+ rw.sd=c(
+ x1.0=1,x2.0=1
+ ),
+ Np=1000,
+ cooling.factor=1,
+ max.fail=10
+ )
+
+print(logLik(pfilter(ou2,Np=2000)),digits=4)
+print(logLik(pfilter(ou2,params=p,Np=2000)),digits=4)
+print(logLik(pfilter(ou2,params=coef(fit),Np=2000)),digits=4)
+print(coef(fit,ics))
+print(coef(ou2,ics))
+print(p-coef(ou2))
+print(coef(fit)-p)
+print(coef(fit)-coef(ou2))
Added: pkg/tests/ou2-icfit.Rout.save
===================================================================
--- pkg/tests/ou2-icfit.Rout.save (rev 0)
+++ pkg/tests/ou2-icfit.Rout.save 2011-09-05 15:37:23 UTC (rev 558)
@@ -0,0 +1,104 @@
+
+R version 2.13.1 (2011-07-08)
+Copyright (C) 2011 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+Platform: x86_64-unknown-linux-gnu (64-bit)
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
+>
+> set.seed(921625222L)
+>
+> data(ou2)
+>
+> ics <- c("x1.0","x2.0") # names of the initial condition parameters
+>
+> p <- coef(ou2)
+> p[ics] <- c(0,0)
+>
+> fit <- mif(
++ ou2,
++ start=p,
++ Nmif=2,
++ ic.lag=1000,
++ var.factor=4,
++ ivps=ics,
++ rw.sd=c(
++ x1.0=1,x2.0=1
++ ),
++ Np=1000,
++ cooling.factor=1,
++ max.fail=10
++ )
+Warning message:
+mif warning: 'ic.lag' = 1000 > 100 = length(time('object')) is nonsensical. Setting 'ic.lag' = 100.
+>
+> fit <- mif(
++ ou2,
++ start=p,
++ Nmif=2,
++ ic.lag=10,
++ var.factor=4,
++ ivps=ics,
++ rw.sd=c(
++ x1.0=1,x2.0=1
++ ),
++ Np=1000,
++ cooling.factor=1,
++ max.fail=10
++ )
+Warning message:
+mif warning: only IVPs are to be estimated, yet 'ic.lag' = 10 < 100 = length(time('object')), so unnecessary work is to be done.
+>
+> fit <- mif(
++ window(ou2,end=10),
++ start=p,
++ Nmif=300,
++ ic.lag=10,
++ var.factor=4,
++ ivps=ics,
++ rw.sd=c(
++ x1.0=1,x2.0=1
++ ),
++ Np=1000,
++ cooling.factor=0.99,
++ max.fail=10
++ )
+>
+> print(logLik(pfilter(ou2,Np=2000)),digits=4)
+[1] -479.4
+> print(logLik(pfilter(ou2,params=p,Np=2000)),digits=4)
+[1] -481.3
+> print(logLik(pfilter(ou2,params=coef(fit),Np=2000)),digits=4)
+[1] -480.3
+> print(coef(fit,ics))
+ x1.0 x2.0
+-2.991659 2.880117
+> print(p-coef(ou2))
+alpha.1 alpha.2 alpha.3 alpha.4 sigma.1 sigma.2 sigma.3 tau x1.0 x2.0
+ 0 0 0 0 0 0 0 0 3 -4
+> print(coef(fit)-p)
+ alpha.1 alpha.2 alpha.3 alpha.4 sigma.1 sigma.2 sigma.3 tau
+ 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
+ x1.0 x2.0
+-2.991659 2.880117
+> print(coef(fit)-coef(ou2))
+ alpha.1 alpha.2 alpha.3 alpha.4 sigma.1 sigma.2
+ 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000
+ sigma.3 tau x1.0 x2.0
+ 0.000000000 0.000000000 0.008340701 -1.119883047
+>
Deleted: pkg/tests/sir-icfit.R
===================================================================
--- pkg/tests/sir-icfit.R 2011-09-03 13:54:05 UTC (rev 557)
+++ pkg/tests/sir-icfit.R 2011-09-05 15:37:23 UTC (rev 558)
@@ -1,153 +0,0 @@
-library(pomp)
-
-set.seed(343435488L)
-
-pdf(file="sir-icfit.pdf")
-
-data(euler.sir)
-po <- window(euler.sir,end=0.2)
-guess <- coef(po)
-ics <- c("S.0","I.0","R.0")
-guess[ics[-3]] <- guess[ics[-3]]+c(0.2,-0.2)
-
-plist <- list(
- probe.marginal("reports",ref=obs(po),order=3,diff=1,transform=sqrt),
- probe.acf("reports",lags=c(0,1,2,3,4,5),transform=sqrt),
- median=probe.median("reports")
- )
-
-summary(pm.true <- probe(po,probes=plist,nsim=100,seed=1066L))
-
-summary(pm.guess <- probe(po,params=guess,probes=plist,nsim=100,seed=1066L))
-
-pm.fit <- probe.match(
- po,
- start=guess,
- probes=plist,
- est=ics[-1],
- method="Nelder-Mead",
- trace=3,
- reltol=1e-5,
- parscale=c(0.1,0.1),
- nsim=100,
- seed=1066L
- )
-
-summary(pm.fit)
-
-comp.table <- cbind(true=exp(coef(po,ics)),guess=exp(guess[ics]),fit=exp(coef(pm.fit,ics)))
-comp.table <- apply(comp.table,2,function(x)x/sum(x))
-comp.table <- rbind(
- comp.table,
- synth.loglik=c(
- summary(pm.true)$synth.loglik,
- summary(pm.guess)$synth.loglik,
- summary(pm.fit)$synth.loglik
- )
- )
-comp.table
-
-x <- sapply(
- list(true=pm.true,guess=pm.guess,fit=pm.fit),
- function (x) trajectory(x,times=time(x),t0=timezero(x))["cases",1,]
- )
-
-plot(range(time(po)),range(c(states(po,"cases"),x)),bty='l',xlab="time",ylab="cases",type='n')
-points(time(po),states(po,"cases"))
-matlines(time(po),x,lty=1,col=c("red","blue","green"))
-legend("topright",lty=c(NA,1,1,1),pch=c(1,NA,NA,NA),bty='n',col=c("black","red","blue","green"),legend=c("actual",colnames(x)))
-
-summary(tm.true <- traj.match(po,eval.only=TRUE))
-
-summary(tm.guess <- traj.match(po,start=guess,eval.only=TRUE))
-
-tm.fit <- traj.match(
- po,
- start=guess,
- est=ics[-1],
- method="sannbox",
- maxit=300,
- trace=2,
- parscale=c(0.1,0.1)
- )
-
-tm.fit <- traj.match(
- tm.fit,
- est=ics[-1],
- method="Nelder-Mead",
- trace=3,
- reltol=1e-8,
- parscale=c(0.1,0.1)
- )
-
-summary(tm.fit)
-
-comp.table <- cbind(true=exp(coef(po,ics)),guess=exp(guess[ics]),fit=exp(coef(tm.fit,ics)))
-comp.table <- apply(comp.table,2,function(x)x/sum(x))
-comp.table <- rbind(
- comp.table,
- loglik=sapply(list(tm.true,tm.guess,tm.fit),logLik)
- )
-comp.table
-
-x <- sapply(
- list(true=tm.true,guess=tm.guess,fit=tm.fit),
- function (x) trajectory(x,times=time(x),t0=timezero(x))["cases",1,]
- )
-
-plot(range(time(po)),range(c(states(po,"cases"),x)),bty='l',xlab="time",ylab="cases",type='n')
-points(time(po),states(po,"cases"))
-matlines(time(po),x,lty=1,col=c("red","blue","green"))
-legend("topright",lty=c(NA,1,1,1),pch=c(1,NA,NA,NA),bty='n',col=c("black","red","blue","green"),legend=c("actual",colnames(x)))
-
-### now try an initial condition fitting approach based on particle filtering
-est <- ics[-1]
-np <- 10000 # number of particles to use
-pp <- array(coef(po),dim=c(length(coef(po)),np),dimnames=list(names(coef(po)),NULL))
-## generate an array of guesses
-guesses <- sobolDesign(lower=guess[est]-0.5,upper=guess[est]+0.5,nseq=np)
-nd <- length(time(po))
-
-## fit the initial conditions using repeated filtering on the initial window of the data
-
-for (j in seq_len(3)) {
- for (k in est) {
- pp[k,] <- guesses[[k]]
- }
- for (k in seq_len(5)) {
- pf <- pfilter(po,params=pp,save.params=TRUE)
- pp <- pf at saved.params[[nd]]
- }
- guesses <- sobolDesign(
- lower=apply(pp[est,],1,min),
- upper=apply(pp[est,],1,max),
- nseq=np
- )
-}
-
-pf.fit <- po
-coef(pf.fit,ics) <- log(apply(apply(exp(pp[ics,]),2,function(x)x/sum(x)),1,mean))
-pf.true <- pfilter(po,Np=2000)
-pf.guess <- pfilter(po,params=guess,Np=2000,max.fail=100)
-pf.fit <- pfilter(pf.fit,Np=2000)
-
-comp.table <- cbind(true=exp(coef(po,ics)),guess=exp(guess[ics]),fit=exp(coef(pf.fit,ics)))
-comp.table <- apply(comp.table,2,function(x)x/sum(x))
-comp.table <- rbind(
- comp.table,
- loglik=sapply(list(pf.true,pf.guess,pf.fit),logLik)
- )
-comp.table
-
-x <- sapply(
- list(true=pf.true,guess=pf.guess,fit=pf.fit),
- function (x) trajectory(x,times=time(x),t0=timezero(x))["cases",1,]
- )
-
-plot(range(time(po)),range(c(states(po,"cases"),x)),bty='l',xlab="time",ylab="cases",type='n')
-points(time(po),states(po,"cases"))
-matlines(time(po),x,lty=1,col=c("red","blue","green"))
-legend("topright",lty=c(NA,1,1,1),pch=c(1,NA,NA,NA),bty='n',col=c("black","red","blue","green"),legend=c("actual",colnames(x)))
-
-dev.off()
-
Deleted: pkg/tests/sir-icfit.Rout.save
===================================================================
--- pkg/tests/sir-icfit.Rout.save 2011-09-03 13:54:05 UTC (rev 557)
+++ pkg/tests/sir-icfit.Rout.save 2011-09-05 15:37:23 UTC (rev 558)
@@ -1,716 +0,0 @@
-
-R Under development (unstable) (2011-07-13 r56369)
-Copyright (C) 2011 The R Foundation for Statistical Computing
-ISBN 3-900051-07-0
-Platform: x86_64-unknown-linux-gnu (64-bit)
-
-R is free software and comes with ABSOLUTELY NO WARRANTY.
-You are welcome to redistribute it under certain conditions.
-Type 'license()' or 'licence()' for distribution details.
-
-R is a collaborative project with many contributors.
-Type 'contributors()' for more information and
-'citation()' on how to cite R or R packages in publications.
-
-Type 'demo()' for some demos, 'help()' for on-line help, or
-'help.start()' for an HTML browser interface to help.
-Type 'q()' to quit R.
-
-> library(pomp)
-Loading required package: mvtnorm
-Loading required package: subplex
-Loading required package: deSolve
->
-> set.seed(343435488L)
->
-> pdf(file="sir-icfit.pdf")
->
-> data(euler.sir)
-> po <- window(euler.sir,end=0.2)
-> guess <- coef(po)
-> ics <- c("S.0","I.0","R.0")
-> guess[ics[-3]] <- guess[ics[-3]]+c(0.2,-0.2)
->
-> plist <- list(
-+ probe.marginal("reports",ref=obs(po),order=3,diff=1,transform=sqrt),
-+ probe.acf("reports",lags=c(0,1,2,3,4,5),transform=sqrt),
-+ median=probe.median("reports")
-+ )
->
-> summary(pm.true <- probe(po,probes=plist,nsim=100,seed=1066L))
-$coef
- gamma mu iota nbasis degree
- 3.258097e+00 -3.912023e+00 -4.605170e+00 3.000000e+00 3.000000e+00
- period beta1 beta2 beta3 beta.sd
- 1.000000e+00 7.090077e+00 7.495542e+00 6.396930e+00 -6.907755e+00
- pop rho S.0 I.0 R.0
- 2.100000e+06 -5.108256e-01 -3.831980e+00 -6.907755e+00 -2.292750e-02
-
-$nsim
-[1] 100
-
-$quantiles
- marg.1 marg.2 marg.3 acf.0.reports acf.1.reports
- 0.96 0.39 0.06 0.45 0.47
-acf.2.reports acf.3.reports acf.4.reports acf.5.reports median
- 0.31 0.20 0.89 0.79 0.59
-
-$pvals
- marg.1 marg.2 marg.3 acf.0.reports acf.1.reports
- 0.0990099 0.7920792 0.1386139 0.9108911 0.9504950
-acf.2.reports acf.3.reports acf.4.reports acf.5.reports median
- 0.6336634 0.4158416 0.2376238 0.4356436 0.8118812
-
-$synth.loglik
-[1] -3.327182
-
->
-> summary(pm.guess <- probe(po,params=guess,probes=plist,nsim=100,seed=1066L))
-$coef
- gamma mu iota nbasis degree
- 3.258097e+00 -3.912023e+00 -4.605170e+00 3.000000e+00 3.000000e+00
- period beta1 beta2 beta3 beta.sd
- 1.000000e+00 7.090077e+00 7.495542e+00 6.396930e+00 -6.907755e+00
- pop rho S.0 I.0 R.0
- 2.100000e+06 -5.108256e-01 -3.631980e+00 -7.107755e+00 -2.292750e-02
-
-$nsim
-[1] 100
-
-$quantiles
- marg.1 marg.2 marg.3 acf.0.reports acf.1.reports
- 0.91 0.28 0.07 0.00 0.00
-acf.2.reports acf.3.reports acf.4.reports acf.5.reports median
- 0.00 0.00 1.00 1.00 0.00
-
-$pvals
- marg.1 marg.2 marg.3 acf.0.reports acf.1.reports
- 0.19801980 0.57425743 0.15841584 0.01980198 0.01980198
-acf.2.reports acf.3.reports acf.4.reports acf.5.reports median
- 0.01980198 0.01980198 0.01980198 0.01980198 0.01980198
-
-$synth.loglik
-[1] -63.04768
-
->
-> pm.fit <- probe.match(
-+ po,
-+ start=guess,
-+ probes=plist,
-+ est=ics[-1],
-+ method="Nelder-Mead",
-+ trace=3,
-+ reltol=1e-5,
-+ parscale=c(0.1,0.1),
-+ nsim=100,
-+ seed=1066L
-+ )
- Nelder-Mead direct search function minimizer
-function value for initial parameters = 63.047679
- Scaled convergence tolerance is 0.000630477
-Stepsize computed as 7.107755
-BUILD 3 4929.990658 63.047679
-HI-REDUCTION 5 601.469297 63.047679
-HI-REDUCTION 7 496.928487 10.499482
-HI-REDUCTION 9 71.488070 10.499482
-HI-REDUCTION 11 63.047679 10.499482
-LO-REDUCTION 13 42.035345 3.573073
-HI-REDUCTION 15 14.564709 3.573073
-HI-REDUCTION 17 10.499482 3.573073
-HI-REDUCTION 19 7.478882 3.573073
-REFLECTION 21 5.096958 2.628314
-HI-REDUCTION 23 3.573073 2.628314
-SHRINK 27 3.617480 2.109975
-LO-REDUCTION 29 2.628314 2.109975
-HI-REDUCTION 31 2.482146 2.109975
-HI-REDUCTION 33 2.203390 1.540292
-SHRINK 37 3.121442 1.540292
-LO-REDUCTION 39 2.324264 1.540292
-SHRINK 43 2.420458 1.540292
-SHRINK 47 4.610749 1.540292
-LO-REDUCTION 49 2.444944 1.540292
-LO-REDUCTION 51 2.439068 1.540292
-HI-REDUCTION 53 2.012399 1.540292
-SHRINK 57 3.233547 1.540292
-LO-REDUCTION 59 1.972827 1.540292
-SHRINK 63 3.404917 1.540292
-LO-REDUCTION 65 2.424849 1.540292
-SHRINK 69 2.424849 1.540292
-HI-REDUCTION 71 2.393003 1.540292
-SHRINK 75 2.393003 1.540292
-SHRINK 79 3.279083 1.540292
-Exiting from Nelder Mead minimizer
- 81 function evaluations used
->
-> summary(pm.fit)
-$coef
- gamma mu iota nbasis degree
- 3.258097e+00 -3.912023e+00 -4.605170e+00 3.000000e+00 3.000000e+00
- period beta1 beta2 beta3 beta.sd
- 1.000000e+00 7.090077e+00 7.495542e+00 6.396930e+00 -6.907755e+00
- pop rho S.0 I.0 R.0
- 2.100000e+06 -5.108256e-01 -3.631980e+00 -6.691678e+00 1.777793e-01
-
-$nsim
-[1] 100
-
-$quantiles
- marg.1 marg.2 marg.3 acf.0.reports acf.1.reports
- 0.89 0.42 0.09 0.52 0.46
-acf.2.reports acf.3.reports acf.4.reports acf.5.reports median
- 0.35 0.23 0.91 0.75 0.43
-
-$pvals
- marg.1 marg.2 marg.3 acf.0.reports acf.1.reports
- 0.2376238 0.8514851 0.1980198 0.9702970 0.9306931
-acf.2.reports acf.3.reports acf.4.reports acf.5.reports median
- 0.7128713 0.4752475 0.1980198 0.5148515 0.8712871
-
-$synth.loglik
-[1] -1.540292
-
-$est
-[1] "I.0" "R.0"
-
-$weights
-[1] 1
-
-$value
-[1] 1.540292
-
-$eval
-[1] 81 NA
-
-$convergence
-[1] 0
-
->
-> comp.table <- cbind(true=exp(coef(po,ics)),guess=exp(guess[ics]),fit=exp(coef(pm.fit,ics)))
-> comp.table <- apply(comp.table,2,function(x)x/sum(x))
-> comp.table <- rbind(
-+ comp.table,
-+ synth.loglik=c(
-+ summary(pm.true)$synth.loglik,
-+ summary(pm.guess)$synth.loglik,
-+ summary(pm.fit)$synth.loglik
-+ )
-+ )
-> comp.table
- true guess fit
-S.0 0.02166667 0.026342137 0.021651353
-I.0 0.00100000 0.000814969 0.001015489
-R.0 0.97733333 0.972842894 0.977333157
-synth.loglik -3.32718185 -63.047679163 -1.540292335
->
-> x <- sapply(
-+ list(true=pm.true,guess=pm.guess,fit=pm.fit),
-+ function (x) trajectory(x,times=time(x),t0=timezero(x))["cases",1,]
-+ )
->
-> plot(range(time(po)),range(c(states(po,"cases"),x)),bty='l',xlab="time",ylab="cases",type='n')
-> points(time(po),states(po,"cases"))
-> matlines(time(po),x,lty=1,col=c("red","blue","green"))
-> legend("topright",lty=c(NA,1,1,1),pch=c(1,NA,NA,NA),bty='n',col=c("black","red","blue","green"),legend=c("actual",colnames(x)))
->
-> summary(tm.true <- traj.match(po,eval.only=TRUE))
-$params
- gamma mu iota nbasis degree
- 3.258097e+00 -3.912023e+00 -4.605170e+00 3.000000e+00 3.000000e+00
- period beta1 beta2 beta3 beta.sd
- 1.000000e+00 7.090077e+00 7.495542e+00 6.396930e+00 -6.907755e+00
- pop rho S.0 I.0 R.0
- 2.100000e+06 -5.108256e-01 -3.831980e+00 -6.907755e+00 -2.292750e-02
-
-$loglik
-[1] -49.09745
-
-$eval
-[1] 1 0
-
-$convergence
-[1] NA
-
-$msg
-[1] "no optimization performed"
-
->
-> summary(tm.guess <- traj.match(po,start=guess,eval.only=TRUE))
-$params
- gamma mu iota nbasis degree
- 3.258097e+00 -3.912023e+00 -4.605170e+00 3.000000e+00 3.000000e+00
- period beta1 beta2 beta3 beta.sd
- 1.000000e+00 7.090077e+00 7.495542e+00 6.396930e+00 -6.907755e+00
- pop rho S.0 I.0 R.0
- 2.100000e+06 -5.108256e-01 -3.631980e+00 -7.107755e+00 -2.292750e-02
-
-$loglik
-[1] -1768.807
-
-$eval
-[1] 1 0
-
-$convergence
-[1] NA
-
-$msg
-[1] "no optimization performed"
-
->
-> tm.fit <- traj.match(
-+ po,
-+ start=guess,
-+ est=ics[-1],
-+ method="sannbox",
-+ maxit=300,
-+ trace=2,
-+ parscale=c(0.1,0.1)
-+ )
-initial evaluation: 1768.807
-iter 1 val= 1768.807 , accept= FALSE
-iter 2 val= 1768.807 , accept= FALSE
-iter 3 val= 1768.807 , accept= FALSE
-iter 4 val= 1768.807 , accept= FALSE
-iter 5 val= 1768.807 , accept= FALSE
-iter 6 val= 1768.807 , accept= FALSE
-iter 7 val= 1768.807 , accept= FALSE
-iter 8 val= 1768.807 , accept= FALSE
-iter 9 val= 1687.885 , accept= TRUE
-iter 10 val= 1687.885 , accept= FALSE
-iter 11 val= 1548.83 , accept= TRUE
-iter 12 val= 1548.83 , accept= FALSE
-iter 13 val= 1548.83 , accept= FALSE
-iter 14 val= 1548.83 , accept= FALSE
-iter 15 val= 1548.83 , accept= FALSE
-iter 16 val= 1548.83 , accept= FALSE
-iter 17 val= 1548.83 , accept= FALSE
-iter 18 val= 1548.83 , accept= FALSE
-iter 19 val= 1548.83 , accept= FALSE
-iter 20 val= 1548.83 , accept= FALSE
-iter 21 val= 1548.83 , accept= FALSE
-iter 22 val= 1548.83 , accept= FALSE
-iter 23 val= 1548.83 , accept= FALSE
-iter 24 val= 1548.83 , accept= FALSE
-iter 25 val= 1548.83 , accept= FALSE
-iter 26 val= 1548.83 , accept= FALSE
-iter 27 val= 1389.144 , accept= TRUE
-iter 28 val= 1389.144 , accept= FALSE
-iter 29 val= 1375.29 , accept= TRUE
-iter 30 val= 976.1314 , accept= TRUE
-iter 31 val= 976.1314 , accept= FALSE
-iter 32 val= 976.1314 , accept= FALSE
-iter 33 val= 976.1314 , accept= FALSE
-iter 34 val= 976.1314 , accept= FALSE
-iter 35 val= 976.1314 , accept= FALSE
-iter 36 val= 976.1314 , accept= FALSE
-iter 37 val= 976.1314 , accept= FALSE
-iter 38 val= 976.1314 , accept= FALSE
-iter 39 val= 976.1314 , accept= FALSE
-iter 40 val= 976.1314 , accept= FALSE
-iter 41 val= 976.1314 , accept= FALSE
-iter 42 val= 976.1314 , accept= FALSE
-iter 43 val= 976.1314 , accept= FALSE
-iter 44 val= 976.1314 , accept= FALSE
-iter 45 val= 976.1314 , accept= FALSE
-iter 46 val= 976.1314 , accept= FALSE
-iter 47 val= 976.1314 , accept= FALSE
-iter 48 val= 976.1314 , accept= FALSE
-iter 49 val= 976.1314 , accept= FALSE
-iter 50 val= 976.1314 , accept= FALSE
-iter 51 val= 976.1314 , accept= FALSE
-iter 52 val= 976.1314 , accept= FALSE
-iter 53 val= 976.1314 , accept= FALSE
-iter 54 val= 976.1314 , accept= FALSE
-iter 55 val= 518.0907 , accept= TRUE
-iter 56 val= 331.9071 , accept= TRUE
-iter 57 val= 331.9071 , accept= FALSE
-iter 58 val= 263.2471 , accept= TRUE
-iter 59 val= 263.2471 , accept= FALSE
-iter 60 val= 263.2471 , accept= FALSE
-iter 61 val= 263.2471 , accept= FALSE
-iter 62 val= 263.2471 , accept= FALSE
-iter 63 val= 263.2471 , accept= FALSE
-iter 64 val= 263.2471 , accept= FALSE
-iter 65 val= 263.2471 , accept= FALSE
-iter 66 val= 263.2471 , accept= FALSE
-iter 67 val= 263.2471 , accept= FALSE
-iter 68 val= 263.2471 , accept= FALSE
-iter 69 val= 263.2471 , accept= FALSE
-iter 70 val= 263.2471 , accept= FALSE
-iter 71 val= 263.2471 , accept= FALSE
-iter 72 val= 263.2471 , accept= FALSE
-iter 73 val= 263.2471 , accept= FALSE
-iter 74 val= 263.2471 , accept= FALSE
-iter 75 val= 263.2471 , accept= FALSE
-iter 76 val= 263.2471 , accept= FALSE
-iter 77 val= 263.2471 , accept= FALSE
-iter 78 val= 263.2471 , accept= FALSE
-iter 79 val= 263.2471 , accept= FALSE
-iter 80 val= 263.2471 , accept= FALSE
-iter 81 val= 263.2471 , accept= FALSE
-iter 82 val= 263.2471 , accept= FALSE
-iter 83 val= 263.2471 , accept= FALSE
-iter 84 val= 263.2471 , accept= FALSE
-iter 85 val= 248.6137 , accept= TRUE
-iter 86 val= 248.6137 , accept= FALSE
-iter 87 val= 248.6137 , accept= FALSE
-iter 88 val= 248.6137 , accept= FALSE
-iter 89 val= 248.6137 , accept= FALSE
-iter 90 val= 248.6137 , accept= FALSE
-iter 91 val= 248.6137 , accept= FALSE
-iter 92 val= 248.6137 , accept= FALSE
-iter 93 val= 248.6137 , accept= FALSE
-iter 94 val= 248.6137 , accept= FALSE
-iter 95 val= 248.6137 , accept= FALSE
-iter 96 val= 240.8026 , accept= TRUE
-iter 97 val= 240.8026 , accept= FALSE
-iter 98 val= 240.8026 , accept= FALSE
-iter 99 val= 240.8026 , accept= FALSE
-iter 100 val= 240.8026 , accept= FALSE
-iter 101 val= 240.8026 , accept= FALSE
-iter 102 val= 240.8026 , accept= FALSE
-iter 103 val= 240.8026 , accept= FALSE
-iter 104 val= 240.8026 , accept= FALSE
-iter 105 val= 240.8026 , accept= FALSE
-iter 106 val= 240.8026 , accept= FALSE
-iter 107 val= 240.8026 , accept= FALSE
-iter 108 val= 240.8026 , accept= FALSE
-iter 109 val= 240.8026 , accept= FALSE
-iter 110 val= 240.8026 , accept= FALSE
-iter 111 val= 240.8026 , accept= FALSE
-iter 112 val= 240.8026 , accept= FALSE
-iter 113 val= 240.8026 , accept= FALSE
-iter 114 val= 240.8026 , accept= FALSE
-iter 115 val= 240.8026 , accept= FALSE
-iter 116 val= 240.8026 , accept= FALSE
-iter 117 val= 240.8026 , accept= FALSE
-iter 118 val= 240.8026 , accept= FALSE
-iter 119 val= 240.8026 , accept= FALSE
-iter 120 val= 240.8026 , accept= FALSE
-iter 121 val= 240.8026 , accept= FALSE
-iter 122 val= 240.8026 , accept= FALSE
-iter 123 val= 240.8026 , accept= FALSE
-iter 124 val= 240.8026 , accept= FALSE
-iter 125 val= 240.8026 , accept= FALSE
-iter 126 val= 240.8026 , accept= FALSE
-iter 127 val= 240.8026 , accept= FALSE
-iter 128 val= 240.8026 , accept= FALSE
-iter 129 val= 240.8026 , accept= FALSE
-iter 130 val= 240.8026 , accept= FALSE
-iter 131 val= 240.8026 , accept= FALSE
-iter 132 val= 240.8026 , accept= FALSE
-iter 133 val= 240.8026 , accept= FALSE
-iter 134 val= 240.8026 , accept= FALSE
-iter 135 val= 240.8026 , accept= FALSE
-iter 136 val= 240.8026 , accept= FALSE
-iter 137 val= 240.8026 , accept= FALSE
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/pomp -r 558
More information about the pomp-commits
mailing list