[Pomp-commits] r568 - in pkg: . R inst inst/doc man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Dec 13 21:25:28 CET 2011
Author: kingaa
Date: 2011-12-13 21:25:25 +0100 (Tue, 13 Dec 2011)
New Revision: 568
Modified:
pkg/DESCRIPTION
pkg/R/bsmc.R
pkg/inst/ChangeLog
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.pdf
pkg/man/bsmc.Rd
pkg/tests/ou2-bsmc.R
pkg/tests/ou2-bsmc.Rout.save
Log:
- fix some inconsistencies with the arguments of 'bsmc'
- switch 'bsmc' to using 'systematic_resampling' instead of 'sample.int'
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2011-11-22 21:49:11 UTC (rev 567)
+++ pkg/DESCRIPTION 2011-12-13 20:25:25 UTC (rev 568)
@@ -2,7 +2,7 @@
Type: Package
Title: Statistical inference for partially observed Markov processes
Version: 0.39-5
-Date: 2011-11-22
+Date: 2011-12-14
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/bsmc.R
===================================================================
--- pkg/R/bsmc.R 2011-11-22 21:49:11 UTC (rev 567)
+++ pkg/R/bsmc.R 2011-12-13 20:25:25 UTC (rev 568)
@@ -4,12 +4,9 @@
##
## params = the initial particles for the parameter values
## these are drawn from the prior distribution for the parameters
-## if a parameter is being held fixed, it is given as a row of NA's
## est = names of parameters to estimate. Other parameters are not updated.
-## Np = number of particles
-## discount = delta, introduced in section 3.3 in Liu & West
+## smooth = parameter 'h' from AGM
## ntries = number of samplesto draw from x_t+1|x(k)_t to estimate mean of mu(k)_t+1 as in sect 2.2 Liu & West
-## ntimes = number of timesteps in observation vector
## lower = lower bounds on prior
## upper = upper bounds on prior
@@ -18,7 +15,7 @@
setMethod(
"bsmc",
"pomp",
- function (object, params, est,
+ function (object, params, Np, est,
smooth = 0.1,
ntries = 1,
tol = 1e-17,
@@ -30,7 +27,8 @@
if (missing(seed)) seed <- NULL
if (!is.null(seed)) {
- if (!exists(".Random.seed",where=.GlobalEnv)) { # need to initialize the RNG
+ if (!exists(".Random.seed",where=.GlobalEnv)) {
+ ## need to initialize the RNG
runif(1)
}
save.seed <- get(".Random.seed",pos=.GlobalEnv)
@@ -47,7 +45,10 @@
}
}
- Np <- NCOL(params)
+ if (missing(Np)) Np <- NCOL(params)
+ else if (is.matrix(params)&&(Np!=ncol(params)))
+ stop(error.prefix,sQuote("Np")," cannot be other than ",sQuote("ncol(params)"),call.=FALSE)
+
ntimes <- length(time(object))
if (is.null(dim(params))) {
params <- matrix(
@@ -64,6 +65,9 @@
npars <- nrow(params)
paramnames <- rownames(params)
+
+ if (missing(est))
+ est <- paramnames[apply(params,1,function(x)diff(range(x))>0)]
estind <- match(est,paramnames)
npars.est <- length(estind)
@@ -138,25 +142,23 @@
)
}
- for (j in seq_len(Np) ) {
- x.tmp <- matrix(data=x[,j],nrow=nvars,ncol=ntries)
- rownames(x.tmp) <- statenames
- p.tmp <- matrix(data=params[,j],nrow=npars,ncol=ntries)
- rownames(p.tmp) <- paramnames
-
- ## update mean of states at time nt as per L&W AGM (1)
- tries <- rprocess(
- object,
- xstart=x.tmp,
- times=times[c(nt,nt+1)],
- params=p.tmp
- )
-
- mu[,j,1] <- apply(tries[,,2,drop=FALSE],1,mean)
-
- ## shrink parameters towards mean as per Liu & West eq (3.3) and L&W AGM (1)
- m[,j] <- shrink*params[,j] + (1-shrink)*params.mean
- }
+ X <- matrix(data=x,nrow=nvars,ncol=Np*ntries)
+ rownames(X) <- statenames
+ P <- matrix(data=params,nrow=npars,ncol=Np*ntries)
+ rownames(P) <- paramnames
+ ## update mean of states at time nt as per L&W AGM (1)
+ tries <- rprocess(
+ object,
+ xstart=X,
+ times=times[c(nt,nt+1)],
+ params=P
+ )[,,2,drop=FALSE]
+ dim(tries) <- c(nvars,Np,ntries,1)
+ mu <- apply(tries,c(1,2,4),mean)
+ rownames(mu) <- statenames
+ ## shrink parameters towards mean as per Liu & West eq (3.3) and L&W AGM (1)
+ m <- shrink*params+(1-shrink)*params.mean
+
## evaluate probability of obervation given mean value of parameters and states (used in L&W AGM (5) below)
g <- dmeasure(
object,
@@ -166,6 +168,7 @@
params=m
)
## sample indices -- From L&W AGM (2)
+## k <- .Call(systematic_resampling,g)
k <- sample.int(n=Np,size=Np,replace=TRUE,prob=g)
params <- params[,k]
m <- m[,k]
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2011-11-22 21:49:11 UTC (rev 567)
+++ pkg/inst/ChangeLog 2011-12-13 20:25:25 UTC (rev 568)
@@ -1,3 +1,14 @@
+2011-11-22 kingaa
+
+ * [r567] R/bsmc.R, inst/NEWS: - update NEWS file
+ * [r566] DESCRIPTION, R/bsmc.R, inst/TODO: - fix bug in 'bsmc' that
+ arises when the state space is 1-D (add drop=FALSE flag to []
+ operation)
+
+2011-11-12 kingaa
+
+ * [r565] DESCRIPTION, inst/ChangeLog, inst/NEWS: - version 0.39-4
+
2011-11-10 kingaa
* [r564] DESCRIPTION, R/sannbox.R, tests/ou2-probe.Rout.save: -
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/man/bsmc.Rd
===================================================================
--- pkg/man/bsmc.Rd 2011-11-22 21:49:11 UTC (rev 567)
+++ pkg/man/bsmc.Rd 2011-12-13 20:25:25 UTC (rev 568)
@@ -8,8 +8,8 @@
\code{bsmc} gives draws from the posterior.
}
\usage{
-\S4method{bsmc}{pomp}(object, params, est, smooth = 0.1, ntries = 1,
- tol = 1e-17, lower = -Inf, upper = Inf, seed = NULL,
+\S4method{bsmc}{pomp}(object, params, Np, est, smooth = 0.1,
+ ntries = 1, tol = 1e-17, lower = -Inf, upper = Inf, seed = NULL,
verbose = getOption("verbose"), max.fail = 0, \dots)
}
\arguments{
@@ -17,15 +17,20 @@
An object of class \code{pomp} or inheriting class \code{pomp}.
}
\item{params}{
- A \code{npars} x \code{np} matrix containing the parameters corresponding to the initial state values in \code{xstart}.
- The matrix should be Np columns long, where Np is the number of particles. The values for each row should be Np draws from the prior distribution for the parameter.
+ A \code{npars} x \code{Np} matrix (with rownames) containing the parameters corresponding to the initial state values in \code{xstart}.
+ \code{Np} is the number of particles, i.e., each row should contain \code{Np} draws from the prior distribution for that parameter.
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, which is equivalent to a degenerate prior.
- If some parameter values are to be held fixed in the fitting (i.e. params.fixed = TRUE), then those elements of params (rows if params is given as a matrix), must be set to NA.
}
+ \item{Np}{
+ If \code{params} is specified as a named vector, \code{Np} specifies the number of particles to use.
+ If \code{params} is specified as a matrix, \code{Np} should not be specified;
+ it is taken to be the number of columns of \code{params}.
+ }
\item{est}{
Names of the rows of \code{params} that are to be estimated.
No updates will be made to the other parameters.
+ If \code{est} is not specified, all parameters for which there is variation in \code{params} will be estimated.
}
\item{smooth}{
Kernel density smoothing parameters.
@@ -37,7 +42,7 @@
Number of draws from \code{rprocess} per particle used to estimate the expected value of the state process at time \code{t+1} given the state and parameters at time \code{t}.
}
\item{tol}{
- Particles with log likelihood below \code{tol} are considered to be "lost".
+ 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)}.
}
Modified: pkg/tests/ou2-bsmc.R
===================================================================
--- pkg/tests/ou2-bsmc.R 2011-11-22 21:49:11 UTC (rev 567)
+++ pkg/tests/ou2-bsmc.R 2011-12-13 20:25:25 UTC (rev 568)
@@ -49,3 +49,8 @@
)
print(min(smc$eff.sample.size))
+print(smc$loglik)
+
+smc <- bsmc(ou2,ntries=5,Np=5000,smooth=0.1,est=estnames,seed=455485L)
+print(smc$eff.sample.size)
+print(smc$loglik)
Modified: pkg/tests/ou2-bsmc.Rout.save
===================================================================
--- pkg/tests/ou2-bsmc.Rout.save 2011-11-22 21:49:11 UTC (rev 567)
+++ pkg/tests/ou2-bsmc.Rout.save 2011-12-13 20:25:25 UTC (rev 568)
@@ -1,5 +1,5 @@
-R version 2.13.1 (2011-07-08)
+R Under development (unstable) (2011-12-13 r57882)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-linux-gnu (64-bit)
@@ -62,7 +62,7 @@
> post <- smc$post
>
> print(etime <- toc-tic)
-Time difference of 38.74774 secs
+Time difference of 3.258152 secs
>
> print(
+ cbind(
@@ -74,8 +74,8 @@
+ )
prior.mean posterior.mean truth 2.5% 50% 97.5%
alpha.1 0.8000000 0.8000000 0.8 0.8000000 0.8000000 0.8000000
-alpha.2 -0.4999287 -0.5043393 -0.5 -0.5112463 -0.5065709 -0.4767260
-alpha.3 0.2996065 0.3270662 0.3 0.3146781 0.3319924 0.3384011
+alpha.2 -0.4999287 -0.4853833 -0.5 -0.5353593 -0.4881214 -0.4521997
+alpha.3 0.2996065 0.3186368 0.3 0.2933113 0.3120256 0.3476582
alpha.4 0.9000000 0.9000000 0.9 0.9000000 0.9000000 0.9000000
sigma.1 3.0000000 3.0000000 3.0 3.0000000 3.0000000 3.0000000
sigma.2 -0.5000000 -0.5000000 -0.5 -0.5000000 -0.5000000 -0.5000000
@@ -85,5 +85,14 @@
x2.0 4.0000000 4.0000000 4.0 4.0000000 4.0000000 4.0000000
>
> print(min(smc$eff.sample.size))
-[1] 6.94596
+[1] 40.66176
+> print(smc$loglik)
+[1] -30.90868
>
+> smc <- bsmc(ou2,ntries=5,Np=5000,smooth=0.1,est=estnames,seed=455485L)
+> print(smc$eff.sample.size)
+ [1] 334.84129 29.59904 21.74851 135.57314 134.38589 13.53892 74.06161
+ [8] 19.83331 30.75126 79.16116
+> print(smc$loglik)
+[1] -30.27596
+>
More information about the pomp-commits
mailing list