[Pomp-commits] r362 - pkg/tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 5 18:27:11 CEST 2010
Author: kingaa
Date: 2010-10-05 18:27:11 +0200 (Tue, 05 Oct 2010)
New Revision: 362
Added:
pkg/tests/ou2-pmcmc.Rout.save
pkg/tests/skeleton.R
pkg/tests/skeleton.Rout.save
Log:
- add in results from ou2-pmcmc test
- add another set of tests in 'skeleton.R'
Added: pkg/tests/ou2-pmcmc.Rout.save
===================================================================
--- pkg/tests/ou2-pmcmc.Rout.save (rev 0)
+++ pkg/tests/ou2-pmcmc.Rout.save 2010-10-05 16:27:11 UTC (rev 362)
@@ -0,0 +1,57 @@
+
+R version 2.11.1 (2010-05-31)
+Copyright (C) 2010 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> library(pomp)
+>
+> data(ou2)
+>
+> dprior.ou2 <- function (params, hyperparams, ..., log) {
++ f <- sum(dunif(params,min=hyperparams$min,max=hyperparams$max,log=TRUE))
++ if (log) f else exp(f)
++ }
+>
+> pdf(file="ou2-pmcmc.pdf")
+>
+> f1 <- pmcmc(
++ ou2,
++ start=coef(ou2),
++ Nmcmc=20,
++ dprior=dprior.ou2,
++ hyperparams=list(min=coef(ou2)-1,max=coef(ou2)+1),
++ rw.sd=c(alpha.2=0.01,alpha.3=0.01),
++ Np=100,
++ max.fail=100,
++ verbose=FALSE
++ )
+> f1 <- continue(f1,Nmcmc=20,max.fail=100)
+> plot(f1)
+>
+> if (FALSE) {
++ f2 <- pmcmc(
++ f1,Nmcmc=1000,Np=500,max.fail=100,
++ verbose=FALSE
++ )
++ plot(f2)
++ runs <- rle(conv.rec(f2)[,'loglik'])$lengths
++ plot(runs)
++ acf(conv.rec(f2)[,c("alpha.2","alpha.3")])
++ }
+>
+> dev.off()
+null device
+ 1
+>
Added: pkg/tests/skeleton.R
===================================================================
--- pkg/tests/skeleton.R (rev 0)
+++ pkg/tests/skeleton.R 2010-10-05 16:27:11 UTC (rev 362)
@@ -0,0 +1,77 @@
+library(pomp)
+
+data(ricker)
+
+x <- array(
+ data=states(ricker),
+ dim=c(2,3,52),
+ dimnames=list(rownames(states(ricker)),NULL,NULL)
+ )
+p <- array(
+ data=coef(ricker),
+ dim=c(5,3),
+ dimnames=list(names(coef(ricker)),NULL)
+ )
+p["log.r",] <- c(1,2,4)
+f <- skeleton(ricker,x=x,params=p,t=time(ricker,t0=T))
+plot(x[1,,],f[1,,],type='n')
+points(x[1,1,],f[1,1,],col='red')
+points(x[1,2,],f[1,2,],col='blue')
+points(x[1,3,],f[1,3,],col='green')
+
+## non-autonomous case
+
+data(euler.sir)
+x <- array(
+ data=states(euler.sir)[,1],
+ dim=c(5,3,100),
+ dimnames=list(rownames(states(euler.sir)),NULL,NULL)
+ )
+p <- array(
+ data=coef(euler.sir),
+ dim=c(15,3),
+ dimnames=list(names(coef(euler.sir)),NULL)
+ )
+p["beta2",1:2] <- c(3,5) ## try different values of one of the seasonality parameters
+## compute the skeleton at each point
+f <- skeleton(euler.sir,x=x,params=p,t=seq(0,1,length=100))
+## verify that the skeleton varies with time
+matplot(seq(0,1,length=100),t(f[1,,]),type='l',lty=1)
+
+x <- trajectory(euler.sir) ## use deSolve to compute a deterministic trajectory
+f <- skeleton(euler.sir,x=x,params=p[,3,drop=F],t=time(euler.sir)) ## evaluate skeleton at each point along it
+
+fit <- smooth.spline(time(euler.sir),x["S",,])## fit a spline to the trajectory
+plot(predict(fit,deriv=1)$y,f["S",,])## compare spline estimate with computed vectorfield
+abline(a=0,b=1)
+
+fit <- smooth.spline(time(euler.sir),x["I",,])
+plot(predict(fit,deriv=1)$y,f["I",,])
+abline(a=0,b=1)
+
+pomp.skeleton <- function(times,y,p,more) {
+# Turns a skeleton function from a 'pomp' object into the right hand
+# side of and ODE for use in CollocInfer
+
+ x <- array(y,dim=c(nrow(y),1,ncol(y)),dimnames=list(rownames(y),NULL,NULL))
+ params <- array(data=p,dim=c(length(p),1),dimnames=list(names(p),NULL))
+
+ skeleton(more$pomp.obj,x=x,params=params,t=times)
+}
+
+x <- array(
+ data=states(ricker),
+ dim=c(2,3,51),
+ dimnames=list(rownames(states(ricker)),NULL,NULL)
+ )
+p <- array(
+ data=coef(ricker),
+ dim=c(5,3),
+ dimnames=list(names(coef(ricker)),NULL)
+ )
+p["log.r",]<- c(1,2,4)
+
+f <- skeleton(ricker,x=x,params=p,t=time(ricker))
+
+pomp.skeleton(time(ricker),x[,1,],p[,1],list(pomp.obj=ricker))
+
Added: pkg/tests/skeleton.Rout.save
===================================================================
--- pkg/tests/skeleton.Rout.save (rev 0)
+++ pkg/tests/skeleton.Rout.save 2010-10-05 16:27:11 UTC (rev 362)
@@ -0,0 +1,406 @@
+
+R version 2.11.1 (2010-05-31)
+Copyright (C) 2010 The R Foundation for Statistical Computing
+ISBN 3-900051-07-0
+
+R is free software and comes with ABSOLUTELY NO WARRANTY.
+You are welcome to redistribute it under certain conditions.
+Type 'license()' or 'licence()' for distribution details.
+
+R is a collaborative project with many contributors.
+Type 'contributors()' for more information and
+'citation()' on how to cite R or R packages in publications.
+
+Type 'demo()' for some demos, 'help()' for on-line help, or
+'help.start()' for an HTML browser interface to help.
+Type 'q()' to quit R.
+
+> library(pomp)
+>
+> data(ricker)
+>
+> x <- array(
++ data=states(ricker),
++ dim=c(2,3,52),
++ dimnames=list(rownames(states(ricker)),NULL,NULL)
++ )
+> p <- array(
++ data=coef(ricker),
++ dim=c(5,3),
++ dimnames=list(names(coef(ricker)),NULL)
++ )
+> p["log.r",] <- c(1,2,4)
+> f <- skeleton(ricker,x=x,params=p,t=time(ricker,t0=T))
+> plot(x[1,,],f[1,,],type='n')
+> points(x[1,1,],f[1,1,],col='red')
+> points(x[1,2,],f[1,2,],col='blue')
+> points(x[1,3,],f[1,3,],col='green')
+>
+> ## non-autonomous case
+>
+> data(euler.sir)
+> x <- array(
++ data=states(euler.sir)[,1],
++ dim=c(5,3,100),
++ dimnames=list(rownames(states(euler.sir)),NULL,NULL)
++ )
+> p <- array(
++ data=coef(euler.sir),
++ dim=c(15,3),
++ dimnames=list(names(coef(euler.sir)),NULL)
++ )
+> p["beta2",1:2] <- c(3,5) ## try different values of one of the seasonality parameters
+> ## compute the skeleton at each point
+> f <- skeleton(euler.sir,x=x,params=p,t=seq(0,1,length=100))
+> ## verify that the skeleton varies with time
+> matplot(seq(0,1,length=100),t(f[1,,]),type='l',lty=1)
+>
+> x <- trajectory(euler.sir) ## use deSolve to compute a deterministic trajectory
+Warning message:
+The default behavior of 'trajectory' has changed.
+See the documentation ("pomp?trajectory") for details
+and set the 'times' and 't0' arguments appropriately to compensate.
+
+> f <- skeleton(euler.sir,x=x,params=p[,3,drop=F],t=time(euler.sir)) ## evaluate skeleton at each point along it
+>
+> fit <- smooth.spline(time(euler.sir),x["S",,])## fit a spline to the trajectory
+> plot(predict(fit,deriv=1)$y,f["S",,])## compare spline estimate with computed vectorfield
+> abline(a=0,b=1)
+>
+> fit <- smooth.spline(time(euler.sir),x["I",,])
+> plot(predict(fit,deriv=1)$y,f["I",,])
+> abline(a=0,b=1)
+>
+> pomp.skeleton <- function(times,y,p,more) {
++ # Turns a skeleton function from a 'pomp' object into the right hand
++ # side of and ODE for use in CollocInfer
++
++ x <- array(y,dim=c(nrow(y),1,ncol(y)),dimnames=list(rownames(y),NULL,NULL))
++ params <- array(data=p,dim=c(length(p),1),dimnames=list(names(p),NULL))
++
++ skeleton(more$pomp.obj,x=x,params=params,t=times)
++ }
+>
+> x <- array(
++ data=states(ricker),
++ dim=c(2,3,51),
++ dimnames=list(rownames(states(ricker)),NULL,NULL)
++ )
+> p <- array(
++ data=coef(ricker),
++ dim=c(5,3),
++ dimnames=list(names(coef(ricker)),NULL)
++ )
+> p["log.r",]<- c(1,2,4)
+>
+> f <- skeleton(ricker,x=x,params=p,t=time(ricker))
+>
+> pomp.skeleton(time(ricker),x[,1,],p[,1],list(pomp.obj=ricker))
+, , 1
+
+ [,1]
+N 0.01735127
+e 0.00000000
+
+, , 2
+
+ [,1]
+N 0.1253023
+e 0.0000000
+
+, , 3
+
+ [,1]
+N 0.0001788570
+e 0.0000000000
+
+, , 4
+
+ [,1]
+N 0.06773162
+e 0.00000000
+
+, , 5
+
+ [,1]
+N 9.646812e-05
+e 0.000000e+00
+
+, , 6
+
+ [,1]
+N 0.4430618
+e 0.0000000
+
+, , 7
+
+ [,1]
+N 0.0009287406
+e 0.0000000000
+
+, , 8
+
+ [,1]
+N 1.117575e-05
+e 0.000000e+00
+
+, , 9
+
+ [,1]
+N 0.6218961
+e 0.0000000
+
+, , 10
+
+ [,1]
+N 0.7720283
+e 0.0000000
+
+, , 11
+
+ [,1]
+N 0.2345794
+e 0.0000000
+
+, , 12
+
+ [,1]
+N 0.01825815
+e 0.00000000
+
+, , 13
+
+ [,1]
+N 0.007791341
+e 0.000000000
+
+, , 14
+
+ [,1]
+N 0.9084504
+e 0.0000000
+
+, , 15
+
+ [,1]
+N 0.7690181
+e 0.0000000
+
+, , 16
+
+ [,1]
+N 0.2664813
+e 0.0000000
+
+, , 17
+
+ [,1]
+N 2.337531e-10
+e 0.000000e+00
+
+, , 18
+
+ [,1]
+N 0.01735127
+e 0.00000000
+
+, , 19
+
+ [,1]
+N 0.1253023
+e 0.0000000
+
+, , 20
+
+ [,1]
+N 0.0001788570
+e 0.0000000000
+
+, , 21
+
+ [,1]
+N 0.06773162
+e 0.00000000
+
+, , 22
+
+ [,1]
+N 9.646812e-05
+e 0.000000e+00
+
+, , 23
+
+ [,1]
+N 0.4430618
+e 0.0000000
+
+, , 24
+
+ [,1]
+N 0.0009287406
+e 0.0000000000
+
+, , 25
+
+ [,1]
+N 1.117575e-05
+e 0.000000e+00
+
+, , 26
+
+ [,1]
+N 0.6218961
+e 0.0000000
+
+, , 27
+
+ [,1]
+N 0.7720283
+e 0.0000000
+
+, , 28
+
+ [,1]
+N 0.2345794
+e 0.0000000
+
+, , 29
+
+ [,1]
+N 0.01825815
+e 0.00000000
+
+, , 30
+
+ [,1]
+N 0.007791341
+e 0.000000000
+
+, , 31
+
+ [,1]
+N 0.9084504
+e 0.0000000
+
+, , 32
+
+ [,1]
+N 0.7690181
+e 0.0000000
+
+, , 33
+
+ [,1]
+N 0.2664813
+e 0.0000000
+
+, , 34
+
+ [,1]
+N 2.337531e-10
+e 0.000000e+00
+
+, , 35
+
+ [,1]
+N 0.01735127
+e 0.00000000
+
+, , 36
+
+ [,1]
+N 0.1253023
+e 0.0000000
+
+, , 37
+
+ [,1]
+N 0.0001788570
+e 0.0000000000
+
+, , 38
+
+ [,1]
+N 0.06773162
+e 0.00000000
+
+, , 39
+
+ [,1]
+N 9.646812e-05
+e 0.000000e+00
+
+, , 40
+
+ [,1]
+N 0.4430618
+e 0.0000000
+
+, , 41
+
+ [,1]
+N 0.0009287406
+e 0.0000000000
+
+, , 42
+
+ [,1]
+N 1.117575e-05
+e 0.000000e+00
+
+, , 43
+
+ [,1]
+N 0.6218961
+e 0.0000000
+
+, , 44
+
+ [,1]
+N 0.7720283
+e 0.0000000
+
+, , 45
+
+ [,1]
+N 0.2345794
+e 0.0000000
+
+, , 46
+
+ [,1]
+N 0.01825815
+e 0.00000000
+
+, , 47
+
+ [,1]
+N 0.007791341
+e 0.000000000
+
+, , 48
+
+ [,1]
+N 0.9084504
+e 0.0000000
+
+, , 49
+
+ [,1]
+N 0.7690181
+e 0.0000000
+
+, , 50
+
+ [,1]
+N 0.2664813
+e 0.0000000
+
+, , 51
+
+ [,1]
+N 2.337531e-10
+e 0.000000e+00
+
+>
+>
More information about the pomp-commits
mailing list