[Pomp-commits] r593 - in pkg: R inst man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jan 12 15:56:13 CET 2012
Author: kingaa
Date: 2012-01-12 15:56:13 +0100 (Thu, 12 Jan 2012)
New Revision: 593
Modified:
pkg/R/bsmc.R
pkg/inst/NEWS
pkg/man/bsmc.Rd
pkg/tests/ou2-bsmc.R
pkg/tests/ou2-bsmc.Rout.save
Log:
- Pierre Jacob's contributions to 'bsmc'
Modified: pkg/R/bsmc.R
===================================================================
--- pkg/R/bsmc.R 2012-01-12 13:54:24 UTC (rev 592)
+++ pkg/R/bsmc.R 2012-01-12 14:56:13 UTC (rev 593)
@@ -2,11 +2,12 @@
##
## in annotation L&W AGM == Liu & West "A General Algorithm"
##
-## params = the initial particles for the parameter values
-## these are drawn from the prior distribution for the parameters
-## est = names of parameters to estimate. Other parameters are not updated.
+## params = the initial particles for the parameter values;
+## these should be drawn from the prior distribution for the parameters
+## est = names of parameters to estimate; other parameters are not updated.
## 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
+## 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
## lower = lower bounds on prior
## upper = upper bounds on prior
@@ -27,10 +28,8 @@
if (missing(seed)) seed <- NULL
if (!is.null(seed)) {
- if (!exists(".Random.seed",where=.GlobalEnv)) {
- ## need to initialize the RNG
- runif(1)
- }
+ if (!exists(".Random.seed",where=.GlobalEnv))
+ runif(1) ## need to initialize the RNG
save.seed <- get(".Random.seed",pos=.GlobalEnv)
set.seed(seed)
}
@@ -116,7 +115,7 @@
times <- time(object,t0=TRUE)
x <- xstart
- loglik <- rep(NA,ntimes)
+ evidence <- rep(NA,ntimes)
eff.sample.size <- rep(NA,ntimes)
nfail <- 0
@@ -142,7 +141,7 @@
)
}
- X <- matrix(data=x,nrow=nvars,ncol=Np*ntries)
+ X <- matrix(data=x,nrow=nvars,ncol=Np*ntries)
rownames(X) <- statenames
P <- matrix(data=params,nrow=npars,ncol=Np*ntries)
rownames(P) <- paramnames
@@ -167,6 +166,7 @@
times=times[nt+1],
params=m
)
+ storeForEvidence1 <- log(sum(g))
## sample indices -- From L&W AGM (2)
## k <- .Call(systematic_resampling,g)
k <- sample.int(n=Np,size=Np,replace=TRUE,prob=g)
@@ -207,7 +207,9 @@
params=params
)
## evaluate weights as per L&W AGM (5)
- weights <- numer/g
+
+ weights <- numer/g
+ storeForEvidence2 <- log(mean(weights))
## apply box constraints as per the priors
for (j in seq_len(Np)) {
@@ -248,12 +250,12 @@
nfail <- nfail+1
if (nfail > max.fail)
stop(error.prefix,"too many filtering failures",call.=FALSE)
- loglik[nt] <- log(tol) # worst log-likelihood
+ evidence[nt] <- log(tol) # worst log-likelihood
weights <- rep(1/Np,Np)
eff.sample.size[nt] <- 0
} else { # not all particles are lost
## compute log-likelihood
- loglik[nt] <- log(mean(weights))
+ evidence[nt] <- storeForEvidence1+storeForEvidence2
weights[failures] <- 0
weights <- weights/sum(weights)
## compute effective sample-size
@@ -283,11 +285,11 @@
post=params,
prior=prior,
eff.sample.size=eff.sample.size,
- cond.loglik=loglik,
smooth=smooth,
seed=seed,
nfail=nfail,
- loglik=sum(loglik),
+ cond.log.evidence=evidence,
+ log.evidence=sum(evidence),
weights=weights
)
}
Modified: pkg/inst/NEWS
===================================================================
--- pkg/inst/NEWS 2012-01-12 13:54:24 UTC (rev 592)
+++ pkg/inst/NEWS 2012-01-12 14:56:13 UTC (rev 593)
@@ -1,4 +1,9 @@
NEWS
+0.40-3
+ o The 'bsmc' method has been improved, due to the contributions of Pierre Jacob.
+ Specifically, 'bsmc' no longer reports a log-likelihood (which it never really computed anyway) but a log-evidence.
+ The latter computation was supplied by Pierre Jacob.
+
0.40-2
o A bug to do with computation of the number of steps needed in discrete-time simulation and trajectory computations has been fixed.
This bug was introduced in version 0.40-1.
Modified: pkg/man/bsmc.Rd
===================================================================
--- pkg/man/bsmc.Rd 2012-01-12 13:54:24 UTC (rev 592)
+++ pkg/man/bsmc.Rd 2012-01-12 14:56:13 UTC (rev 593)
@@ -78,9 +78,6 @@
\item{eff.sample.size}{
A vector containing the effective number of particles at each time point.
}
- \item{cond.loglik}{
- A vector containing the conditional log likelihoods at each time point.
- }
\item{smooth}{
The smoothing parameter used (see above).
}
@@ -92,9 +89,12 @@
\item{nfail}{
The number of filtering failures encountered.
}
- \item{loglik}{
- The estimated log-likelihood.
+ \item{cond.log.evidence}{
+ A vector containing the conditional log evidence scores at each time point.
}
+ \item{log.evidence}{
+ The estimated log evidence.
+ }
\item{weights}{
The resampling weights for each particle.
}
Modified: pkg/tests/ou2-bsmc.R
===================================================================
--- pkg/tests/ou2-bsmc.R 2012-01-12 13:54:24 UTC (rev 592)
+++ pkg/tests/ou2-bsmc.R 2012-01-12 14:56:13 UTC (rev 593)
@@ -49,8 +49,8 @@
)
print(min(smc$eff.sample.size))
-print(smc$loglik)
+print(smc$log.evidence)
smc <- bsmc(ou2,ntries=5,Np=5000,smooth=0.1,est=estnames,seed=455485L)
print(smc$eff.sample.size)
-print(smc$loglik)
+print(smc$log.evidence)
Modified: pkg/tests/ou2-bsmc.Rout.save
===================================================================
--- pkg/tests/ou2-bsmc.Rout.save 2012-01-12 13:54:24 UTC (rev 592)
+++ pkg/tests/ou2-bsmc.Rout.save 2012-01-12 14:56:13 UTC (rev 593)
@@ -1,5 +1,5 @@
-R Under development (unstable) (2011-12-13 r57882)
+R version 2.14.1 (2011-12-22)
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 3.258152 secs
+Time difference of 3.09274 secs
>
> print(
+ cbind(
@@ -86,13 +86,13 @@
>
> print(min(smc$eff.sample.size))
[1] 40.66176
-> print(smc$loglik)
-[1] -30.90868
+> print(smc$log.evidence)
+[1] 44.55899
>
> 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
+> print(smc$log.evidence)
+[1] 38.46819
>
More information about the pomp-commits
mailing list