[Pomp-commits] r403 - in pkg: . R inst man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Oct 25 12:58:28 CEST 2010
Author: kingaa
Date: 2010-10-25 12:58:28 +0200 (Mon, 25 Oct 2010)
New Revision: 403
Modified:
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/spect-match.R
pkg/R/spect.R
pkg/inst/NEWS
pkg/man/mif.Rd
pkg/man/probe.Rd
pkg/tests/ou2-mif.R
pkg/tests/ou2-mif.Rout.save
pkg/tests/ricker.R
pkg/tests/ricker.Rout.save
Log:
- remove 'seed' slot from 'spect.pomp' objects
- tinker with tests
- 'tests/ricker.R' now contains a MIF vs PM contest
- temporarily remove 'mif.profile.design' from the visible interface
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2010-10-24 15:14:57 UTC (rev 402)
+++ pkg/DESCRIPTION 2010-10-25 10:58:28 UTC (rev 403)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.34-3
-Date: 2010-10-24
+Version: 0.35-1
+Date: 2010-10-25
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>
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-10-24 15:14:57 UTC (rev 402)
+++ pkg/NAMESPACE 2010-10-25 10:58:28 UTC (rev 403)
@@ -65,7 +65,6 @@
bspline.basis,
periodic.bspline.basis,
compare.mif,
- mif.profile.design,
nlf,
traj.match,
probe.mean,
Modified: pkg/R/spect-match.R
===================================================================
--- pkg/R/spect-match.R 2010-10-24 15:14:57 UTC (rev 402)
+++ pkg/R/spect-match.R 2010-10-25 10:58:28 UTC (rev 403)
@@ -85,6 +85,8 @@
verbose = getOption("verbose"),
eval.only = FALSE, fail.value = NA, ...) {
+ obj.fn <- spect.mismatch
+
if (!is(object,"pomp"))
stop(sQuote("object")," must be of class ",sQuote("pomp"))
@@ -156,21 +158,21 @@
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
- )
+ val <- obj.fn(
+ 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")
@@ -178,7 +180,7 @@
if (method == 'subplex') {
opt <- subplex::subplex(
par=guess,
- fn=spect.mismatch,
+ fn=obj.fn,
est=par.index,
object=object,
params=params,
@@ -196,7 +198,7 @@
} else {
opt <- optim(
par=guess,
- fn=spect.mismatch,
+ fn=obj.fn,
est=par.index,
object=object,
params=params,
Modified: pkg/R/spect.R
===================================================================
--- pkg/R/spect.R 2010-10-24 15:14:57 UTC (rev 402)
+++ pkg/R/spect.R 2010-10-25 10:58:28 UTC (rev 403)
@@ -7,7 +7,6 @@
"spect.pomp",
contains="pomp",
representation=representation(
- seed="integer",
kernel.width="numeric",
transform="function",
freq="numeric",
@@ -159,12 +158,6 @@
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,
@@ -209,7 +202,6 @@
new(
"spect.pomp",
object,
- seed=as.integer(seed),
kernel.width=kernel.width,
transform=transform,
detrend=detrend,
@@ -229,11 +221,6 @@
if (missing(vars)) probes <- rownames(object at datspec)
if (missing(kernel.width)) kernel.width <- object at kernel.width
if (missing(nsim)) nsim <- nrow(object at simspec)
- if (is.null(seed)) {
- if (exists('.Random.seed',where=.GlobalEnv)) {
- seed <- get(".Random.seed",pos=.GlobalEnv)
- }
- }
if (missing(transform)) transform <- object at transform
if (missing(detrend)) detrend <- object at detrend
spect(
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2010-10-24 15:14:57 UTC (rev 402)
+++ pkg/inst/NEWS 2010-10-25 10:58:28 UTC (rev 403)
@@ -1,11 +1,19 @@
-NEW IN VERSION 0.34-3:
+NEW IN VERSION 0.35-1:
o Added a facility for computation of "synthetic likelihood" (Wood 2010) using a robust estimate of the covariance. A new slot in objects of class 'probed.pomp' holds this value.
+ o The synthetic likelihood is now the quantity that is maximized in 'probe.match'.
+
o Added the argument 'type' to 'probe.ccf' so that either cross-covariance or cross-correlation can be returned. 'probe.ccf' and 'probe.acf' now have comparable interfaces.
o Changed the way that 'probe.acf' accepts lags. Now an arbitrary set of lags can be selected.
+ o A bug in the new method 'obs' has been fixed.
+
+ o A protection-stack overflow bug in 'probe' has been fixed.
+
+ o The 'seed' slot in 'probed.pomp' and 'spect.pomp' objects has been removed, along with unwanted behavior to do with it.
+
NEW IN VERSION 0.34-2:
o Added a method to coerce 'probed.pomp' objects to 'data.frame' objects.
Modified: pkg/man/mif.Rd
===================================================================
--- pkg/man/mif.Rd 2010-10-24 15:14:57 UTC (rev 402)
+++ pkg/man/mif.Rd 2010-10-25 10:58:28 UTC (rev 403)
@@ -8,7 +8,7 @@
\alias{continue}
\alias{continue,mif-method}
\alias{continue-mif}
-\alias{mif.profile.design}
+%%\alias{mif.profile.design}
\title{The MIF algorithm}
\description{
The MIF algorithm for estimating the parameters of a partially-observed Markov process.
@@ -21,8 +21,8 @@
verbose = getOption("verbose"))
\S4method{mif}{mif}(object, Nmif, \dots)
\S4method{continue}{mif}(object, Nmif = 1, \dots)
-mif.profile.design(object, profile, lower, upper, nprof, ivps, rw.sd,
- Np, ic.lag, var.factor, cooling.factor, \dots)
+%%mif.profile.design(object, profile, lower, upper, nprof, ivps, rw.sd,
+%% Np, ic.lag, var.factor, cooling.factor, \dots)
}
\arguments{
\item{object}{
@@ -92,14 +92,14 @@
\item{verbose}{
logical; if TRUE, print progress reports.
}
- \item{lower, upper, nprof}{
- used in specifying the profile calculation.
- See \code{\link{profile.design}} for details.
- }
- \item{profile}{
- named list of numeric vectors.
- The profile calculation will fix parameters at all possible combinations of the parameters in \code{profile}.
- }
+% \item{lower, upper, nprof}{
+% used in specifying the profile calculation.
+% See \code{\link{profile.design}} for details.
+% }
+% \item{profile}{
+% named list of numeric vectors.
+% The profile calculation will fix parameters at all possible combinations of the parameters in \code{profile}.
+% }
\item{\dots}{
Additional arguments that can be used to override the defaults.
}
@@ -152,7 +152,7 @@
\author{Aaron A. King \email{kingaa at umich dot edu}}
\seealso{
\code{\link{mif-class}}, \code{\link{mif-methods}}, \code{\link{pomp}}, \code{\link{pomp-class}}, \code{\link{pfilter}}.
- See the \dQuote{intro_to_pomp} vignette for an example.
+ See the \dQuote{intro_to_pomp} vignette for examples.
}
\keyword{models}
\keyword{ts}
Modified: pkg/man/probe.Rd
===================================================================
--- pkg/man/probe.Rd 2010-10-24 15:14:57 UTC (rev 402)
+++ pkg/man/probe.Rd 2010-10-25 10:58:28 UTC (rev 403)
@@ -14,6 +14,7 @@
\code{probe} applies one or more \dQuote{probes} to time series data and model simulations and compares the results.
It can be used to diagnose goodness of fit and/or as the basis for \dQuote{probe-matching}, a generalized method-of-moments approach to parameter estimation.
\code{probe.match} calls an optimizer to adjust model parameters to do probe-matching, i.e., to minimize the discrepancy between simulated and actual data.
+ This discrepancy is measured using the \dQuote{synthetic likelihood} as defined by Wood (2010).
}
\usage{
probe(object, probes, \dots)
@@ -93,6 +94,10 @@
B. E. Kendall, C. J. Briggs, W. M. Murdoch, P. Turchin, S. P. Ellner, E. McCauley, R. M. Nisbet, S. N. Wood
Why do populations cycle? A synthesis of statistical and mechanistic modeling approaches,
Ecology, 80:1789--1805, 1999.
+
+ S. N. Wood
+ Statistical inference for noisy nonlinear ecological dynamic systems,
+ Nature, 466: 1102--1104, 2010.
}
\author{
Daniel C. Reuman (d.reuman at imperial dot ac dot uk)
Modified: pkg/tests/ou2-mif.R
===================================================================
--- pkg/tests/ou2-mif.R 2010-10-24 15:14:57 UTC (rev 402)
+++ pkg/tests/ou2-mif.R 2010-10-25 10:58:28 UTC (rev 403)
@@ -184,17 +184,17 @@
)
pp <- particles(fit,Np=10,center=coef(fit),sd=abs(0.1*coef(fit)))
-mp <- mif.profile.design(
- ou2,
- profile=list(
- alpha.1=seq(0.5,0.9,by=0.1),
- alpha.4=seq(0.5,0.9,by=0.1),
- sigma.1=1,sigma.2=0,sigma.3=1,
- x1.0=0,x2.0=0
- ),
- lower=c(alpha.2=-1,alpha.3=-1,tau=0.01),
- upper=c(alpha.2=1,alpha.3=1,tau=3),
- nprof=5,
- rw.sd=c(alpha.2=0.2,alpha.3=0.2,tau=0.1),
- Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
- )
+## mp <- mif.profile.design(
+## ou2,
+## profile=list(
+## alpha.1=seq(0.5,0.9,by=0.1),
+## alpha.4=seq(0.5,0.9,by=0.1),
+## sigma.1=1,sigma.2=0,sigma.3=1,
+## x1.0=0,x2.0=0
+## ),
+## lower=c(alpha.2=-1,alpha.3=-1,tau=0.01),
+## upper=c(alpha.2=1,alpha.3=1,tau=3),
+## nprof=5,
+## rw.sd=c(alpha.2=0.2,alpha.3=0.2,tau=0.1),
+## Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
+## )
Modified: pkg/tests/ou2-mif.Rout.save
===================================================================
--- pkg/tests/ou2-mif.Rout.save 2010-10-24 15:14:57 UTC (rev 402)
+++ pkg/tests/ou2-mif.Rout.save 2010-10-25 10:58:28 UTC (rev 403)
@@ -246,18 +246,18 @@
+ )
> pp <- particles(fit,Np=10,center=coef(fit),sd=abs(0.1*coef(fit)))
>
-> mp <- mif.profile.design(
-+ ou2,
-+ profile=list(
-+ alpha.1=seq(0.5,0.9,by=0.1),
-+ alpha.4=seq(0.5,0.9,by=0.1),
-+ sigma.1=1,sigma.2=0,sigma.3=1,
-+ x1.0=0,x2.0=0
-+ ),
-+ lower=c(alpha.2=-1,alpha.3=-1,tau=0.01),
-+ upper=c(alpha.2=1,alpha.3=1,tau=3),
-+ nprof=5,
-+ rw.sd=c(alpha.2=0.2,alpha.3=0.2,tau=0.1),
-+ Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
-+ )
+> ## mp <- mif.profile.design(
+> ## ou2,
+> ## profile=list(
+> ## alpha.1=seq(0.5,0.9,by=0.1),
+> ## alpha.4=seq(0.5,0.9,by=0.1),
+> ## sigma.1=1,sigma.2=0,sigma.3=1,
+> ## x1.0=0,x2.0=0
+> ## ),
+> ## lower=c(alpha.2=-1,alpha.3=-1,tau=0.01),
+> ## upper=c(alpha.2=1,alpha.3=1,tau=3),
+> ## nprof=5,
+> ## rw.sd=c(alpha.2=0.2,alpha.3=0.2,tau=0.1),
+> ## Np=100,cooling.factor=0.95,ic.lag=10,var.factor=1
+> ## )
>
Modified: pkg/tests/ricker.R
===================================================================
--- pkg/tests/ricker.R 2010-10-24 15:14:57 UTC (rev 402)
+++ pkg/tests/ricker.R 2010-10-25 10:58:28 UTC (rev 403)
@@ -2,16 +2,34 @@
data(ricker)
-set.seed(64857673L)
+pdf(file="ricker.pdf")
+tj.1 <- trajectory(ricker)
+plot(time(ricker),tj.1[1,,],type='l')
+tj.2 <- trajectory(ricker,times=c(30:50),t0=0)
+lines(30:50,tj.2[1,,],col='red',lwd=2)
+max(abs(tj.1[,,time(ricker)>=30]-tj.2[,,]))
+
po <- ricker
-coef(po) <- c(log.r=log(2),log.sigma=log(1),log.phi=log(20),N.0=7,e.0=0)
+try(
+ coef(po,"r")
+ )
+coef(po,c("r","phi")) <- c(0,0)
+coef(po,c("log.r","log.phi")) <- c(a=0,b=0)
+coef(po,c("log.r","log.phi")) <- 0
+coef(po) <- c(log.phi=0,log.r=3.5,N.0=10,e.0=0,log.sigma=-Inf)
+coef(po)
+coef(po,"new") <- 3
+plot(simulate(po))
+coef(po)
-pf.tr <- pfilter(ricker,Np=1000,max.fail=50)
-pf.po <- pfilter(po,Np=1000,max.fail=50)
+set.seed(64857673L)
+guess <- ricker
+coef(guess) <- c(log.r=log(20),log.sigma=log(1),log.phi=log(20),N.0=7,e.0=0)
+
mf <- mif(
- po,
+ guess,
Nmif=100,
Np=1000,
cooling.factor=0.99,
@@ -20,32 +38,55 @@
max.fail=50,
rw.sd=c(log.r=0.1,log.sigma=0.1,log.phi=0.1)
)
-plot(mf <- mif(mf,Nmif=500,max.fail=20))
-pf.mf <- pfilter(mf,Np=1000)
+mf <- continue(mf,Nmif=500,max.fail=20)
+plist <- list(
+ probe.marginal("y",ref=obs(ricker),transform=sqrt),
+ probe.acf("y",lags=c(0,1,2,3,4),transform=sqrt),
+ probe.nlar("y",lags=c(1,1,1,2),powers=c(1,2,3,1),transform=sqrt)
+ )
+
+pr.guess <- probe(guess,probes=plist,nsim=1000,seed=1066L)
+pr.tr <- probe(ricker,probes=plist,nsim=1000,seed=1066L)
+
+pm <- probe.match(
+ pr.guess,
+ est=c("log.r","log.sigma","log.phi"),
+ method="Nelder-Mead",
+ maxit=2000,
+ seed=1066L,
+ reltol=1e-8,
+ trace=3
+ )
+
+pf.tr <- pfilter(ricker,Np=1000,max.fail=50,seed=1066L)
+pf.guess <- pfilter(guess,Np=1000,max.fail=50,seed=1066L)
+pf.mf <- pfilter(mf,Np=1000,seed=1066L)
+pf.pm <- pfilter(pm,Np=1000,max.fail=10,seed=1066L)
+
+pr.mf <- probe(mf,nsim=1000,probes=plist,seed=1066L)
+
res <- rbind(
- cbind(truth=coef(ricker),MLE=coef(mf),guess=coef(po)),
- loglik=c(pf.tr$loglik,pf.mf$loglik,pf.po$loglik)
+ cbind(guess=coef(guess),truth=coef(ricker),MLE=coef(mf),PM=coef(pm)),
+ loglik=c(
+ pf.guess$loglik,
+ pf.tr$loglik,
+ pf.mf$loglik,
+ pf.pm$loglik
+ ),
+ synth.loglik=c(
+ summary(pr.guess)$synth.loglik,
+ summary(pr.tr)$synth.loglik,
+ summary(pr.mf)$synth.loglik,
+ summary(pm)$synth.loglik
+ )
)
print(res,digits=3)
-tj.1 <- trajectory(ricker)
-plot(time(ricker),tj.1[1,,],type='l')
-tj.2 <- trajectory(ricker,times=c(30:50),t0=0)
-lines(30:50,tj.2[1,,],col='red',lwd=2)
-max(abs(tj.1[,,time(ricker)>=30]-tj.2[,,]))
+plot(ricker)
+plot(simulate(guess))
+plot(simulate(mf))
+plot(simulate(pm))
-data(ricker)
-po <- ricker
-try(
- coef(po,"r")
- )
-coef(po,c("r","phi")) <- c(0,0)
-coef(po,c("log.r","log.phi")) <- c(a=0,b=0)
-coef(po,c("log.r","log.phi")) <- 0
-coef(po) <- c(log.phi=0,log.r=3.5,N.0=10,e.0=0,log.sigma=-Inf)
-coef(po)
-coef(po,"new") <- 3
-plot(simulate(po))
-coef(po)
+dev.off()
Modified: pkg/tests/ricker.Rout.save
===================================================================
--- pkg/tests/ricker.Rout.save 2010-10-24 15:14:57 UTC (rev 402)
+++ pkg/tests/ricker.Rout.save 2010-10-25 10:58:28 UTC (rev 403)
@@ -19,41 +19,8 @@
>
> data(ricker)
>
-> set.seed(64857673L)
+> pdf(file="ricker.pdf")
>
-> po <- ricker
-> coef(po) <- c(log.r=log(2),log.sigma=log(1),log.phi=log(20),N.0=7,e.0=0)
->
-> pf.tr <- pfilter(ricker,Np=1000,max.fail=50)
-> pf.po <- pfilter(po,Np=1000,max.fail=50)
->
-> mf <- mif(
-+ po,
-+ Nmif=100,
-+ Np=1000,
-+ cooling.factor=0.99,
-+ var.factor=2,
-+ ic.lag=3,
-+ max.fail=50,
-+ rw.sd=c(log.r=0.1,log.sigma=0.1,log.phi=0.1)
-+ )
-> plot(mf <- mif(mf,Nmif=500,max.fail=20))
-> pf.mf <- pfilter(mf,Np=1000)
->
-> res <- rbind(
-+ cbind(truth=coef(ricker),MLE=coef(mf),guess=coef(po)),
-+ loglik=c(pf.tr$loglik,pf.mf$loglik,pf.po$loglik)
-+ )
->
-> print(res,digits=3)
- truth MLE guess
-log.r 3.80 3.81 0.693
-log.sigma -1.20 -1.73 0.000
-log.phi 2.30 2.32 2.996
-N.0 7.00 7.00 7.000
-e.0 0.00 0.00 0.000
-loglik -139.22 -137.14 -268.810
->
> tj.1 <- trajectory(ricker)
Warning message:
The default behavior of 'trajectory' has changed.
@@ -66,7 +33,6 @@
> max(abs(tj.1[,,time(ricker)>=30]-tj.2[,,]))
[1] 0
>
-> data(ricker)
> po <- ricker
> try(
+ coef(po,"r")
@@ -92,3 +58,213 @@
log.phi log.r N.0 e.0 log.sigma new
0.0 3.5 10.0 0.0 -Inf 3.0
>
+> set.seed(64857673L)
+>
+> guess <- ricker
+> coef(guess) <- c(log.r=log(20),log.sigma=log(1),log.phi=log(20),N.0=7,e.0=0)
+>
+> mf <- mif(
++ guess,
++ Nmif=100,
++ Np=1000,
++ cooling.factor=0.99,
++ var.factor=2,
++ ic.lag=3,
++ max.fail=50,
++ rw.sd=c(log.r=0.1,log.sigma=0.1,log.phi=0.1)
++ )
+> mf <- continue(mf,Nmif=500,max.fail=20)
+>
+> plist <- list(
++ probe.marginal("y",ref=obs(ricker),transform=sqrt),
++ probe.acf("y",lags=c(0,1,2,3,4),transform=sqrt),
++ probe.nlar("y",lags=c(1,1,1,2),powers=c(1,2,3,1),transform=sqrt)
++ )
+>
+> pr.guess <- probe(guess,probes=plist,nsim=1000,seed=1066L)
+> pr.tr <- probe(ricker,probes=plist,nsim=1000,seed=1066L)
+>
+> pm <- probe.match(
++ pr.guess,
++ est=c("log.r","log.sigma","log.phi"),
++ method="Nelder-Mead",
++ maxit=2000,
++ seed=1066L,
++ reltol=1e-8,
++ trace=3
++ )
+ Nelder-Mead direct search function minimizer
+function value for initial parameters = 12.281731
+ Scaled convergence tolerance is 1.22817e-07
+Stepsize computed as 0.299573
+BUILD 4 57.293574 12.281731
+LO-REDUCTION 6 50.563448 12.281731
+EXTENSION 8 22.273083 -5.683156
+EXTENSION 10 21.238197 -12.104646
+LO-REDUCTION 12 12.281731 -12.104646
+LO-REDUCTION 14 -5.683156 -12.104646
+HI-REDUCTION 16 -9.299555 -12.104646
+EXTENSION 18 -10.238701 -16.111147
+LO-REDUCTION 20 -10.432496 -16.111147
+LO-REDUCTION 22 -12.104646 -16.111147
+EXTENSION 24 -13.678744 -17.527973
+EXTENSION 26 -14.637732 -18.135507
+LO-REDUCTION 28 -16.111147 -18.406162
+REFLECTION 30 -17.527973 -18.901071
+HI-REDUCTION 32 -18.135507 -19.025728
+HI-REDUCTION 34 -18.406162 -19.025728
+REFLECTION 36 -18.822244 -19.711127
+HI-REDUCTION 38 -18.901071 -19.711127
+HI-REDUCTION 40 -19.025728 -19.711127
+SHRINK 45 -18.987087 -19.858319
+LO-REDUCTION 47 -19.288693 -19.858319
+LO-REDUCTION 49 -19.711127 -19.858319
+LO-REDUCTION 51 -19.779236 -19.858319
+SHRINK 56 -19.098736 -19.858319
+LO-REDUCTION 58 -19.213348 -19.858319
+LO-REDUCTION 60 -19.230757 -19.858319
+SHRINK 65 -19.572673 -19.858319
+HI-REDUCTION 67 -19.622339 -19.858319
+SHRINK 72 -19.417796 -19.883729
+LO-REDUCTION 74 -19.441471 -19.883729
+HI-REDUCTION 76 -19.458238 -19.883729
+LO-REDUCTION 78 -19.622452 -19.883729
+LO-REDUCTION 80 -19.692344 -19.883729
+SHRINK 85 -19.276517 -19.883729
+HI-REDUCTION 87 -19.332435 -19.883729
+LO-REDUCTION 89 -19.391462 -19.883729
+LO-REDUCTION 91 -19.766486 -19.883729
+SHRINK 96 -19.508802 -19.883729
+LO-REDUCTION 98 -19.685618 -19.883729
+SHRINK 103 -19.279074 -19.883729
+SHRINK 108 -19.498156 -20.036526
+SHRINK 113 -19.163003 -20.036526
+LO-REDUCTION 115 -19.685799 -20.036526
+SHRINK 120 -19.437896 -20.036526
+HI-REDUCTION 122 -19.594982 -20.036526
+LO-REDUCTION 124 -19.655126 -20.036526
+HI-REDUCTION 126 -19.693313 -20.036526
+HI-REDUCTION 128 -19.710519 -20.036526
+SHRINK 133 -19.212997 -20.036526
+LO-REDUCTION 135 -19.632972 -20.036526
+SHRINK 140 -19.318310 -20.036526
+LO-REDUCTION 142 -19.429013 -20.036526
+LO-REDUCTION 144 -19.541768 -20.036526
+SHRINK 149 -19.499454 -20.036526
+HI-REDUCTION 151 -19.643647 -20.036526
+HI-REDUCTION 153 -19.665508 -20.036526
+LO-REDUCTION 155 -19.842500 -20.036526
+SHRINK 160 -19.714203 -20.036526
+LO-REDUCTION 162 -19.900941 -20.036526
+SHRINK 167 -19.494052 -20.036526
+LO-REDUCTION 169 -19.557532 -20.036526
+LO-REDUCTION 171 -19.635954 -20.036526
+SHRINK 176 -19.571857 -20.036526
+HI-REDUCTION 178 -19.633935 -20.036526
+LO-REDUCTION 180 -19.773673 -20.036526
+SHRINK 185 -19.580741 -20.036526
+LO-REDUCTION 187 -19.713589 -20.036526
+SHRINK 192 -19.778851 -20.075344
+LO-REDUCTION 194 -19.920306 -20.075344
+SHRINK 199 -19.253544 -20.075344
+LO-REDUCTION 201 -19.563461 -20.075344
+LO-REDUCTION 203 -19.657461 -20.075344
+SHRINK 208 -19.386987 -20.075344
+LO-REDUCTION 210 -19.445799 -20.075344
+LO-REDUCTION 212 -19.595884 -20.075344
+SHRINK 217 -19.582926 -20.151824
+LO-REDUCTION 219 -19.753819 -20.151824
+LO-REDUCTION 221 -19.772159 -20.151824
+LO-REDUCTION 223 -19.782335 -20.151824
+SHRINK 228 -19.596088 -20.151824
+LO-REDUCTION 230 -19.681504 -20.151824
+SHRINK 235 -19.609725 -20.151824
+HI-REDUCTION 237 -19.615594 -20.151824
+HI-REDUCTION 239 -19.651915 -20.151824
+LO-REDUCTION 241 -19.694368 -20.151824
+LO-REDUCTION 243 -19.769449 -20.151824
+HI-REDUCTION 245 -19.782970 -20.151824
+SHRINK 250 -19.592628 -20.151824
+LO-REDUCTION 252 -19.877327 -20.151824
+HI-REDUCTION 254 -19.966120 -20.151824
+SHRINK 259 -19.757861 -20.151824
+LO-REDUCTION 261 -19.978390 -20.151824
+SHRINK 266 -19.866279 -20.151824
+HI-REDUCTION 268 -20.048979 -20.151824
+LO-REDUCTION 270 -20.073212 -20.164817
+HI-REDUCTION 272 -20.082314 -20.265947
+LO-REDUCTION 274 -20.151824 -20.265947
+SHRINK 279 -20.086588 -20.265947
+SHRINK 284 -20.054913 -20.265947
+REFLECTION 286 -20.079536 -20.312392
+LO-REDUCTION 288 -20.216102 -20.312392
+SHRINK 293 -20.153482 -20.312392
+SHRINK 298 -20.097052 -20.312392
+LO-REDUCTION 300 -20.118531 -20.312392
+LO-REDUCTION 302 -20.165018 -20.312392
+LO-REDUCTION 304 -20.176834 -20.312392
+SHRINK 309 -20.130573 -20.312392
+LO-REDUCTION 311 -20.175138 -20.312392
+HI-REDUCTION 313 -20.232404 -20.312392
+SHRINK 318 -20.158381 -20.340133
+LO-REDUCTION 320 -20.250809 -20.340133
+SHRINK 325 -20.236229 -20.340133
+LO-REDUCTION 327 -20.270202 -20.340133
+LO-REDUCTION 329 -20.313522 -20.340133
+SHRINK 334 -20.301190 -20.340133
+LO-REDUCTION 336 -20.320514 -20.340133
+SHRINK 341 -20.309539 -20.340133
+LO-REDUCTION 343 -20.328223 -20.340133
+SHRINK 348 -20.316445 -20.341834
+LO-REDUCTION 350 -20.337619 -20.341834
+HI-REDUCTION 352 -20.340133 -20.341834
+HI-REDUCTION 354 -20.341082 -20.341834
+SHRINK 359 -20.341834 -20.343672
+SHRINK 364 -20.341834 -20.343672
+SHRINK 369 -20.341834 -20.343672
+Polytope size measure not decreased in shrink
+Exiting from Nelder Mead minimizer
+ 371 function evaluations used
+>
+> pf.tr <- pfilter(ricker,Np=1000,max.fail=50,seed=1066L)
+> pf.guess <- pfilter(guess,Np=1000,max.fail=50,seed=1066L)
+> pf.mf <- pfilter(mf,Np=1000,seed=1066L)
+> pf.pm <- pfilter(pm,Np=1000,max.fail=10,seed=1066L)
+>
+> pr.mf <- probe(mf,nsim=1000,probes=plist,seed=1066L)
+>
+> res <- rbind(
++ cbind(guess=coef(guess),truth=coef(ricker),MLE=coef(mf),PM=coef(pm)),
++ loglik=c(
++ pf.guess$loglik,
++ pf.tr$loglik,
++ pf.mf$loglik,
++ pf.pm$loglik
++ ),
++ synth.loglik=c(
++ summary(pr.guess)$synth.loglik,
++ summary(pr.tr)$synth.loglik,
++ summary(pr.mf)$synth.loglik,
++ summary(pm)$synth.loglik
++ )
++ )
+>
+> print(res,digits=3)
+ guess truth MLE PM
+log.r 3.00 3.80 3.81 3.74
+log.sigma 0.00 -1.20 -1.65 -1.12
+log.phi 3.00 2.30 2.32 2.41
+N.0 7.00 7.00 7.00 7.00
+e.0 0.00 0.00 0.00 0.00
+loglik -230.86 -139.55 -137.17 -143.33
+synth.loglik -12.28 17.47 17.59 20.34
+>
+> plot(ricker)
+> plot(simulate(guess))
+> plot(simulate(mf))
+> plot(simulate(pm))
+>
+> dev.off()
+null device
+ 1
+>
More information about the pomp-commits
mailing list