[Pomp-commits] r102 - in pkg: . data inst/doc tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 24 19:34:23 CEST 2009
Author: kingaa
Date: 2009-04-24 19:34:23 +0200 (Fri, 24 Apr 2009)
New Revision: 102
Added:
pkg/data/ou2.R
Removed:
pkg/data/ou2.rda
Modified:
pkg/DESCRIPTION
pkg/inst/doc/compiled_code_in_pomp.Rnw
pkg/tests/ou2-kalman.Rout.save
pkg/tests/ou2-mif.Rout.save
Log:
replace the binary data 'ou2.rda' with R code in 'ou2.R'
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2009-04-21 11:04:16 UTC (rev 101)
+++ pkg/DESCRIPTION 2009-04-24 17:34:23 UTC (rev 102)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.23-2
-Date: 2009-04-17
+Version: 0.23-3
+Date: 2009-04-24
Author: Aaron A. King, Edward L. Ionides, Carles Martinez Breto, Steve Ellner, Bruce Kendall
Maintainer: Aaron A. King <kingaa at umich.edu>
Description: Inference methods for partially-observed Markov processes
Added: pkg/data/ou2.R
===================================================================
--- pkg/data/ou2.R (rev 0)
+++ pkg/data/ou2.R 2009-04-24 17:34:23 UTC (rev 102)
@@ -0,0 +1,77 @@
+require(pomp)
+ou2 <- pomp(
+ times=seq(1,100),
+ data=rbind(
+ y1=rep(0,100),
+ y2=rep(0,100)
+ ),
+ t0=0,
+ rprocess = function (xstart, times, params, paramnames, ...) {
+ nvar <- nrow(xstart)
+ npar <- nrow(params)
+ nrep <- ncol(xstart)
+ ntimes <- length(times)
+ ## get indices of the various parameters in the 'params' matrix
+ ## C uses zero-based indexing!
+ parindex <- match(paramnames,rownames(params))-1
+ array(
+ .C("ou2_adv",
+ X = double(nvar*nrep*ntimes),
+ xstart = as.double(xstart),
+ par = as.double(params),
+ times = as.double(times),
+ n = as.integer(c(nvar,npar,nrep,ntimes)),
+ parindex = as.integer(parindex),
+ DUP = FALSE,
+ NAOK = TRUE,
+ PACKAGE = "pomp"
+ )$X,
+ dim=c(nvar,nrep,ntimes),
+ dimnames=list(rownames(xstart),NULL,NULL)
+ )
+ },
+ dprocess = function (x, times, params, log, ...) {
+ nvar <- nrow(x)
+ npar <- nrow(params)
+ nrep <- ncol(x)
+ ntimes <- length(times)
+ parindex <- match(paramnames,rownames(params))-1
+ array(
+ .C("ou2_pdf",
+ d = double(nrep*(ntimes-1)),
+ X = as.double(x),
+ par = as.double(params),
+ times = as.double(times),
+ n = as.integer(c(nvar,npar,nrep,ntimes)),
+ parindex = as.integer(parindex),
+ give_log=as.integer(log),
+ DUP = FALSE,
+ NAOK = TRUE,
+ PACKAGE = "pomp"
+ )$d,
+ dim=c(nrep,ntimes-1)
+ )
+ },
+ dmeasure = "normal_dmeasure",
+ rmeasure = "normal_rmeasure",
+ obsnames = c("y1","y2"),
+ paramnames = c(
+ "alpha.1","alpha.2","alpha.3","alpha.4",
+ "sigma.1","sigma.2","sigma.3",
+ "tau"
+ ),
+ statenames = c("x1","x2")
+ )
+
+ou2 <- simulate(
+ ou2,
+ params=c(
+ alpha.1=0.9,alpha.2=0,alpha.3=0,alpha.4=0.99,
+ sigma.1=1,sigma.2=0,sigma.3=2,
+ tau=1,x1.0=50,x2.0=-50
+ ),
+ nsim=1,
+ seed=377456545L
+ )
+
+ou2 <- ou2[[1]]
Deleted: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/compiled_code_in_pomp.Rnw
===================================================================
--- pkg/inst/doc/compiled_code_in_pomp.Rnw 2009-04-21 11:04:16 UTC (rev 101)
+++ pkg/inst/doc/compiled_code_in_pomp.Rnw 2009-04-24 17:34:23 UTC (rev 102)
@@ -158,6 +158,7 @@
\begin{figure}
<<fig=T,echo=F>>=
+ data(ou2)
plot(ou2)
@
\caption{One realization of the two-dimensional OU process.}
@@ -165,8 +166,5 @@
\end{figure}
The \code{pomp} object we just created is included in the package: use \code{data(ou2)} to retrieve it.
-<<echo=F>>=
-save(list='ou2',file='ou2.rda')
-@
\end{document}
Modified: pkg/tests/ou2-kalman.Rout.save
===================================================================
--- pkg/tests/ou2-kalman.Rout.save 2009-04-21 11:04:16 UTC (rev 101)
+++ pkg/tests/ou2-kalman.Rout.save 2009-04-24 17:34:23 UTC (rev 102)
@@ -1,6 +1,6 @@
-R version 2.6.1 (2007-11-26)
-Copyright (C) 2007 The R Foundation for Statistical Computing
+R version 2.8.1 (2008-12-22)
+Copyright (C) 2008 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
R is free software and comes with ABSOLUTELY NO WARRANTY.
@@ -16,6 +16,8 @@
Type 'q()' to quit R.
> library(pomp)
+Loading required package: deSolve
+Loading required package: subplex
>
> kalman.filter <- function (y, x0, a, b, sigma, tau) {
+ n <- nrow(y)
@@ -80,7 +82,7 @@
> cat("Kalman filter log likelihood at truth\n")
Kalman filter log likelihood at truth
> print(-kalman(p.truth,ou2,p.truth),digits=4)
-[1] -422.4
+[1] -412
>
> # make a wild guess
> p.guess <- c(alpha.1=0.8,alpha.4=0.9,x1.0=45,x2.0=-60)
@@ -92,7 +94,7 @@
> cat("Kalman filter log likelihood at guess\n")
Kalman filter log likelihood at guess
> print(-kalman(p.guess,ou2,p.truth),digits=4)
-[1] -562.3
+[1] -668.4
>
> # find MLE using Kalman filter starting at the guess
> cat("running Kalman filter estimation\n")
@@ -101,10 +103,10 @@
> kalm.fit1 <- optim(p.guess,kalman,object=ou2,params=p.truth,hessian=T)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 28.87264 secs
+Time difference of 39.45996 secs
> tic <- Sys.time()
> print(-kalm.fit1$value,digits=4)
-[1] -420.7
+[1] -411.1
>
> cat("summary of results\n")
summary of results
@@ -117,8 +119,8 @@
+ digits=4
+ )
truth MLE SE
-alpha.1 0.90 0.9020 0.009790
-alpha.4 0.99 0.9767 0.008197
-x1.0 50.00 50.7439 1.651465
-x2.0 -50.00 -49.3766 2.305126
+alpha.1 0.90 0.9021 0.009832
+alpha.4 0.99 0.9952 0.005294
+x1.0 50.00 50.4747 1.649801
+x2.0 -50.00 -47.9919 2.229328
>
Modified: pkg/tests/ou2-mif.Rout.save
===================================================================
--- pkg/tests/ou2-mif.Rout.save 2009-04-21 11:04:16 UTC (rev 101)
+++ pkg/tests/ou2-mif.Rout.save 2009-04-24 17:34:23 UTC (rev 102)
@@ -32,7 +32,7 @@
> cat("particle filter log likelihood at truth\n")
particle filter log likelihood at truth
> print(fit1.pfilter$loglik,digits=4)
-[1] -423.3
+[1] -412.9
>
> p.truth <- coef(ou2)
> guess <- p.truth
@@ -47,7 +47,7 @@
> cat("particle filter log likelihood at guess\n")
particle filter log likelihood at guess
> print(fit2.pfilter$loglik,digits=4)
-[1] -686.8
+[1] -3913
>
> cat("running MIF\n")
running MIF
@@ -67,11 +67,11 @@
> mif.fit <- continue(mif.fit,Nmif=70,max.fail=100)
> toc <- Sys.time()
> print(toc-tic)
-Time difference of 58.27099 secs
+Time difference of 1.020114 mins
> cat("PF estimated log likelihood at MIF MLE\n")
PF estimated log likelihood at MIF MLE
> print(pfilter(mif.fit)$loglik,digits=4)
-[1] -421.7
+[1] -411.2
>
> cat("coefficients at truth\n")
coefficients at truth
@@ -82,7 +82,7 @@
MIF MLE
> print(coef(mif.fit,c('x1.0','x2.0','alpha.1','alpha.4')),digits=4)
x1.0 x2.0 alpha.1 alpha.4
- 51.2028 -50.0965 0.9032 0.9738
+ 50.8435 -48.7910 0.9035 0.9902
>
> plot(mif.fit)
> compare.mif(mif.fit)
@@ -118,7 +118,7 @@
> fit <- continue(fit,Nmif=40)
> ff <- pfilter(fit,pred.mean=T,filter.mean=T,pred.var=T,max.fail=100)
> print(ff$loglik)
-[1] -477.4053
+[1] -441.1279
> fit <- mif(fit,rw.sd=c(x1.0=5,x2.0=5,alpha.1=0.1,alpha.4=0.1))
> fit <- continue(fit,Nmif=2,ivps=c("x1.0"),pars=c("alpha.1"))
Warning message:
More information about the pomp-commits
mailing list