[Pomp-commits] r392 - in pkg: . R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Oct 20 22:42:58 CEST 2010
Author: kingaa
Date: 2010-10-20 22:42:57 +0200 (Wed, 20 Oct 2010)
New Revision: 392
Modified:
pkg/NAMESPACE
pkg/R/mif.R
pkg/R/probe-match.R
pkg/R/probe.R
pkg/man/mif.Rd
pkg/man/probed-pomp-class.Rd
pkg/tests/ou2-mif.R
Log:
- remove 'seed' slot from 'probed.pomp'-class objects
- add alternative 'neg.synth.loglik' objective function for 'probe-match'
- add new 'mif.profile.design' function for setting up 'mif' profiling computations
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2010-10-20 18:01:15 UTC (rev 391)
+++ pkg/NAMESPACE 2010-10-20 20:42:57 UTC (rev 392)
@@ -65,6 +65,7 @@
bspline.basis,
periodic.bspline.basis,
compare.mif,
+ mif.profile.design,
nlf,
traj.match,
probe.mean,
Modified: pkg/R/mif.R
===================================================================
--- pkg/R/mif.R 2010-10-20 18:01:15 UTC (rev 391)
+++ pkg/R/mif.R 2010-10-20 20:42:57 UTC (rev 392)
@@ -384,3 +384,45 @@
obj
}
)
+
+mif.profile.design <- function (object, profile, lower, upper, nprof, ivps,
+ rw.sd, Np, ic.lag, var.factor, cooling.factor, ...)
+ {
+ if (missing(profile)) profile <- list()
+ if (missing(lower)) lower <- numeric(0)
+ if (missing(upper)) upper <- lower
+ if (length(lower)!=length(upper))
+ stop(sQuote("lower")," and ",sQuote("upper")," must be of the same length")
+ pars <- names(lower)
+ if (missing(ivps)) ivps <- character(0)
+ Np <- as.integer(Np)
+
+ pd <- do.call(profile.design,c(profile,list(lower=lower,upper=upper,nprof=nprof)))
+
+ object <- as(object,"pomp")
+
+ pp <- coef(object)
+ idx <- !(names(pp)%in%names(pd))
+ if (any(idx)) pd <- cbind(pd,as.list(pp[idx]))
+
+ ans <- vector(mode="list",length=nrow(pd))
+ for (k in seq_len(nrow(pd))) {
+ ans[[k]] <- list(
+ mf=mif(
+ object,
+ Nmif=0,
+ start=unlist(pd[k,]),
+ pars=pars,
+ ivps=ivps,
+ rw.sd=rw.sd,
+ Np=Np,
+ ic.lag=ic.lag,
+ var.factor=var.factor,
+ cooling.factor=cooling.factor,
+ ...
+ )
+ )
+ }
+
+ ans
+ }
Modified: pkg/R/probe-match.R
===================================================================
--- pkg/R/probe-match.R 2010-10-20 18:01:15 UTC (rev 391)
+++ pkg/R/probe-match.R 2010-10-20 20:42:57 UTC (rev 392)
@@ -66,6 +66,31 @@
mismatch
}
+neg.synth.loglik <- function (par, est, object, probes, params,
+ nsim = 1, seed = NULL,
+ weights, datval,
+ fail.value = NA) {
+ if (missing(par)) par <- numeric(0)
+ if (missing(est)) est <- integer(0)
+ if (missing(params)) params <- coef(object)
+
+ params[est] <- par
+
+ ## apply probes to model simulations
+ simval <- .Call(
+ apply_probe_sim,
+ object=object,
+ nsim=nsim,
+ params=params,
+ seed=seed,
+ probes=probes,
+ datval=datval
+ )
+
+ ll <- .Call(synth_loglik,simval,datval)
+ -ll
+}
+
probe.match <- function(object, start, est = character(0),
probes, weights,
nsim, seed = NULL,
@@ -73,6 +98,8 @@
verbose = getOption("verbose"),
eval.only = FALSE, fail.value = NA, ...) {
+ obj.fn <- probe.mismatch
+
if (!is(object,"pomp"))
stop(sQuote("object")," must be of class ",sQuote("pomp"))
@@ -107,26 +134,26 @@
datval <- .Call(apply_probe_data,object,probes) # apply probes to data
if (eval.only) {
- val <- probe.mismatch(
- par=guess,
- est=par.index,
- object=object,
- probes=probes,
- params=params,
- nsim=nsim,
- seed=seed,
- weights=weights,
- datval=datval,
- fail.value=fail.value
- )
+ val <- obj.fun(
+ par=guess,
+ est=par.index,
+ object=object,
+ probes=probes,
+ params=params,
+ nsim=nsim,
+ seed=seed,
+ weights=weights,
+ datval=datval,
+ fail.value=fail.value
+ )
conv <- 0
evals <- as.integer(c(1,0))
- msg <- paste(sQuote("probe.mismatch"),"evaluated")
+ msg <- paste("no optimization performed")
} else {
if (method == 'subplex') {
opt <- subplex::subplex(
par=guess,
- fn=probe.mismatch,
+ fn=obj.fn,
est=par.index,
object=object,
probes=probes,
@@ -141,7 +168,7 @@
} else {
opt <- optim(
par=guess,
- fn=probe.mismatch,
+ fn=obj.fn,
est=par.index,
object=object,
probes=probes,
@@ -164,7 +191,13 @@
new(
"probe.matched.pomp",
- probe(object,probes=probes,params=params,nsim=nsim,seed=seed),
+ probe(
+ as(object,"pomp"),
+ probes=probes,
+ params=params,
+ nsim=nsim,
+ seed=seed
+ ),
weights=weights,
fail.value=as.numeric(fail.value),
value=val,
Modified: pkg/R/probe.R
===================================================================
--- pkg/R/probe.R 2010-10-20 18:01:15 UTC (rev 391)
+++ pkg/R/probe.R 2010-10-20 20:42:57 UTC (rev 392)
@@ -3,7 +3,6 @@
contains="pomp",
representation(
probes="list",
- seed="integer",
datvals="numeric",
simvals="array",
quantiles="numeric",
@@ -41,12 +40,6 @@
if (missing(params)) params <- coef(object)
- if (is.null(seed)) {
- if (exists('.Random.seed',where=.GlobalEnv)) {
- seed <- get(".Random.seed",pos=.GlobalEnv)
- }
- }
-
## apply probes to data
datval <- .Call(apply_probe_data,object,probes)
## apply probes to model simulations
@@ -75,12 +68,11 @@
ll <- .Call(synth_loglik,simval,datval)
coef(object) <- params
-
+
new(
"probed.pomp",
object,
probes=probes,
- seed=as.integer(seed),
datvals=datval,
simvals=simval,
quantiles=quants,
@@ -102,12 +94,6 @@
if (missing(params)) params <- coef(object)
- if (is.null(seed)) {
- if (exists('.Random.seed',where=.GlobalEnv)) {
- seed <- get(".Random.seed",pos=.GlobalEnv)
- }
- }
-
if (missing(nsim)) nsim <- nrow(object at simvals)
probe(
Modified: pkg/man/mif.Rd
===================================================================
--- pkg/man/mif.Rd 2010-10-20 18:01:15 UTC (rev 391)
+++ pkg/man/mif.Rd 2010-10-20 20:42:57 UTC (rev 392)
@@ -8,6 +8,7 @@
\alias{continue}
\alias{continue,mif-method}
\alias{continue-mif}
+\alias{mif.profile.design}
\title{The MIF algorithm}
\description{
The MIF algorithm for estimating the parameters of a partially-observed Markov process.
@@ -20,6 +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)
}
\arguments{
\item{object}{
@@ -89,6 +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{\dots}{
Additional arguments that can be used to override the defaults.
}
@@ -104,6 +115,13 @@
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{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.
+ Running \code{mif} or \code{continue} on each of these will maximize the likelihood over the desired parameters, while holding others fixed.
+ Each \code{mif} calculation is independent of the others: the computation is therefore readily parallelized.
+ The list returned by \code{mif.profile.design} will have one element for every row in \code{do.call(expand.grid,profile)}.
+}
\section{Details}{
If \code{particles} is not specified, the default behavior is to draw the particles from a multivariate normal distribution.
\strong{It is the user's responsibility to ensure that, if the optional \code{particles} argument is given, that the \code{particles} function satisfies the following conditions:}
Modified: pkg/man/probed-pomp-class.Rd
===================================================================
--- pkg/man/probed-pomp-class.Rd 2010-10-20 18:01:15 UTC (rev 391)
+++ pkg/man/probed-pomp-class.Rd 2010-10-20 20:42:57 UTC (rev 392)
@@ -18,7 +18,6 @@
A full description of slots in a \code{probed.pomp} or \code{probe.matched.pomp} object follows.
\describe{
\item{probes}{list of the probes applied.}
- \item{seed}{the seed of the RNG used.}
\item{datvals, simvals}{values of each of the probes applied to the real and simulated data, respectively.}
\item{quantiles}{fraction of simulations with probe values less than the value of the probe of the data.}
\item{pvals}{two-sided p-values: fraction of the \code{simvals} that deviate more extremely from the mean of the \code{simvals} than does \code{datavals}.}
Modified: pkg/tests/ou2-mif.R
===================================================================
--- pkg/tests/ou2-mif.R 2010-10-20 18:01:15 UTC (rev 391)
+++ pkg/tests/ou2-mif.R 2010-10-20 20:42:57 UTC (rev 392)
@@ -4,6 +4,9 @@
set.seed(64857673L)
+obs(window(ou2,end=20,start=15))
+obs(window(ou2,end=5),"y1")
+
fit1.pfilter <- pfilter(ou2,Np=1000)
cat("coefficients at `truth'\n")
print(coef(ou2,c('x1.0','x2.0','alpha.2','alpha.3')),digits=4)
@@ -180,3 +183,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
+ )
More information about the pomp-commits
mailing list