[Pomp-commits] r178 - pkg/data
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Nov 16 12:18:56 CET 2009
Author: kingaa
Date: 2009-11-16 12:18:55 +0100 (Mon, 16 Nov 2009)
New Revision: 178
Modified:
pkg/data/euler.sir.R
pkg/data/ou2.R
pkg/data/rw2.R
pkg/data/verhulst.R
Log:
- cosmetology
Modified: pkg/data/euler.sir.R
===================================================================
--- pkg/data/euler.sir.R 2009-11-15 17:26:12 UTC (rev 177)
+++ pkg/data/euler.sir.R 2009-11-16 11:18:55 UTC (rev 178)
@@ -1,56 +1,56 @@
require(pomp)
-euler.sir <- local(
- {
- po <- pomp(
- times=seq(1/52,4,by=1/52),
- data=rbind(measles=numeric(52*4)),
- t0=0,
- tcovar=seq(0,50,by=1/52),
- covar=matrix(
- periodic.bspline.basis(seq(0,50,by=1/52),nbasis=3,period=1,degree=3),
- ncol=3,
- dimnames=list(NULL,paste("seas",1:3,sep=''))
- ),
- delta.t=1/52/20,
- statenames=c("S","I","R","cases","W","B","dW"),
- paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
- covarnames=c("seas1"),
- zeronames=c("cases"),
- comp.names=c("S","I","R"),
- step.fun="sir_euler_simulator",
- rprocess=euler.simulate,
- dens.fun="sir_euler_density",
- dprocess=onestep.density,
- skeleton.vectorfield="sir_ODE",
- rmeasure="binom_rmeasure",
- dmeasure="binom_dmeasure",
- PACKAGE="pomp",
- initializer=function(params, t0, comp.names, ...){
- p <- exp(params)
- snames <- c(
- "S","I","R","cases","W","B",
- "SI","SD","IR","ID","RD","dW"
- )
- fracs <- p[paste(comp.names,"0",sep=".")]
- x0 <- numeric(length(snames))
- names(x0) <- snames
- x0[comp.names] <- round(p['pop']*fracs/sum(fracs))
- x0
- }
- )
-
- coef(po) <- log(
- c(
- gamma=26,mu=0.02,iota=0.01,
- beta1=1200,beta2=1800,beta3=600,
- beta.sd=1e-3,
- pop=2.1e6,
- rho=0.6,
- S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
- )
- )
-
- simulate(po,nsim=1,seed=329348545L)
+local(
+ {
+ po <- pomp(
+ times=seq(1/52,4,by=1/52),
+ data=rbind(measles=numeric(52*4)),
+ t0=0,
+ tcovar=seq(0,50,by=1/52),
+ covar=matrix(
+ periodic.bspline.basis(seq(0,50,by=1/52),nbasis=3,period=1,degree=3),
+ ncol=3,
+ dimnames=list(NULL,paste("seas",1:3,sep=''))
+ ),
+ delta.t=1/52/20,
+ statenames=c("S","I","R","cases","W","B","dW"),
+ paramnames=c("gamma","mu","iota","beta1","beta.sd","pop","rho"),
+ covarnames=c("seas1"),
+ zeronames=c("cases"),
+ comp.names=c("S","I","R"),
+ step.fun="sir_euler_simulator",
+ rprocess=euler.simulate,
+ dens.fun="sir_euler_density",
+ dprocess=onestep.density,
+ skeleton.vectorfield="sir_ODE",
+ rmeasure="binom_rmeasure",
+ dmeasure="binom_dmeasure",
+ PACKAGE="pomp",
+ initializer=function(params, t0, comp.names, ...){
+ p <- exp(params)
+ snames <- c(
+ "S","I","R","cases","W","B",
+ "SI","SD","IR","ID","RD","dW"
+ )
+ fracs <- p[paste(comp.names,"0",sep=".")]
+ x0 <- numeric(length(snames))
+ names(x0) <- snames
+ x0[comp.names] <- round(p['pop']*fracs/sum(fracs))
+ x0
}
)
+
+ coef(po) <- log(
+ c(
+ gamma=26,mu=0.02,iota=0.01,
+ beta1=1200,beta2=1800,beta3=600,
+ beta.sd=1e-3,
+ pop=2.1e6,
+ rho=0.6,
+ S.0=26/1200,I.0=0.001,R.0=1-0.001-26/1200
+ )
+ )
+
+ simulate(po,nsim=1,seed=329348545L)
+ }
+ ) -> euler.sir
Modified: pkg/data/ou2.R
===================================================================
--- pkg/data/ou2.R 2009-11-15 17:26:12 UTC (rev 177)
+++ pkg/data/ou2.R 2009-11-16 11:18:55 UTC (rev 178)
@@ -1,76 +1,76 @@
require(pomp)
-ou2 <- local(
- {
- po <- 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, paramnames, ...) {
- 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",
- paramnames = c(
- "alpha.1","alpha.2","alpha.3","alpha.4",
- "sigma.1","sigma.2","sigma.3",
- "tau"
- ),
- statenames = c("x1","x2")
- )
-
- coef(po) <- 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
- )
-
- simulate(po,nsim=1,seed=377456545L)
- }
- )
+local(
+ {
+ po <- 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, paramnames, ...) {
+ 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",
+ paramnames = c(
+ "alpha.1","alpha.2","alpha.3","alpha.4",
+ "sigma.1","sigma.2","sigma.3",
+ "tau"
+ ),
+ statenames = c("x1","x2")
+ )
+
+ coef(po) <- 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
+ )
+
+ simulate(po,nsim=1,seed=377456545L)
+ }
+ ) -> ou2
Modified: pkg/data/rw2.R
===================================================================
--- pkg/data/rw2.R 2009-11-15 17:26:12 UTC (rev 177)
+++ pkg/data/rw2.R 2009-11-16 11:18:55 UTC (rev 178)
@@ -1,39 +1,39 @@
require(pomp)
-rw2 <- local(
- {
- po <- pomp(
- rprocess = onestep.simulate,
- dprocess = onestep.density,
- step.fun = function(x, t, params, delta.t, ...) {
- c(
- x1=rnorm(n=1,mean=x['x1'],sd=params['s1']*delta.t),
- x2=rnorm(n=1,mean=x['x2'],sd=params['s2']*delta.t)
- )
- },
- dens.fun = function (x1, t1, x2, t2, params, ...) {
- sum(
- dnorm(
- x=x2[c('x1','x2')],
- mean=x1[c('x1','x2')],
- sd=params[c('s1','s2')]*(t2-t1),
- log=TRUE
- ),
- na.rm=TRUE
- )
- },
- measurement.model=list(
- y1 ~ norm(mean=x1,sd=tau),
- y2 ~ norm(mean=x2,sd=tau)
- ),
- times=1:100,
- data=rbind(
- y1=rep(0,100),
- y2=rep(0,100)
- ),
- t0=0
- )
+local(
+ {
+ po <- pomp(
+ rprocess = onestep.simulate,
+ dprocess = onestep.density,
+ step.fun = function(x, t, params, delta.t, ...) {
+ c(
+ x1=rnorm(n=1,mean=x['x1'],sd=params['s1']*delta.t),
+ x2=rnorm(n=1,mean=x['x2'],sd=params['s2']*delta.t)
+ )
+ },
+ dens.fun = function (x1, t1, x2, t2, params, ...) {
+ sum(
+ dnorm(
+ x=x2[c('x1','x2')],
+ mean=x1[c('x1','x2')],
+ sd=params[c('s1','s2')]*(t2-t1),
+ log=TRUE
+ ),
+ na.rm=TRUE
+ )
+ },
+ measurement.model=list(
+ y1 ~ norm(mean=x1,sd=tau),
+ y2 ~ norm(mean=x2,sd=tau)
+ ),
+ times=1:100,
+ data=rbind(
+ y1=rep(0,100),
+ y2=rep(0,100)
+ ),
+ t0=0
+ )
- simulate(po,params=c(x1.0=0,x2.0=0,s1=1,s2=3,tau=1),nsim=1,seed=738377475L)
- }
- )
+ simulate(po,params=c(x1.0=0,x2.0=0,s1=1,s2=3,tau=1),nsim=1,seed=738377475L)
+ }
+ ) -> rw2
Modified: pkg/data/verhulst.R
===================================================================
--- pkg/data/verhulst.R 2009-11-15 17:26:12 UTC (rev 177)
+++ pkg/data/verhulst.R 2009-11-16 11:18:55 UTC (rev 178)
@@ -1,49 +1,54 @@
require(pomp)
-verhulst <- simulate(
- pomp(
- data=rbind(obs=rep(0,1000)),
- times=seq(0.1,by=0.1,length=1000),
- t0=0,
- rprocess=euler.simulate,
- step.fun=function(x,t,params,delta.t,...){
- with(
- as.list(c(x,params)),
- rnorm(
- n=1,
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t)
- )
- )
- },
- dprocess=onestep.density,
- dens.fun=function(x1,x2,t1,t2,params,log,...){
- delta.t <- t2-t1
- with(
- as.list(c(x1,params)),
- dnorm(
- x=x2['n'],
- mean=n+r*n*(1-n/K)*delta.t,
- sd=sigma*n*sqrt(delta.t),
- log=log
- )
- )
- },
- measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
- skeleton.vectorfield=function(x,t,params,...){
- with(
- as.list(c(x,params)),
- r*n*(1-n/K)
- )
- },
- delta.t=0.01
- ),
- params=c(
- n.0=10000,
- K=10000,
- r=0.9,
- sigma=0.4,
- tau=0.1
- ),
- seed=73658676L
- )
+local(
+ {
+ po <- pomp(
+ data=rbind(obs=rep(0,1000)),
+ times=seq(0.1,by=0.1,length=1000),
+ t0=0,
+ rprocess=euler.simulate,
+ step.fun=function(x,t,params,delta.t,...){
+ with(
+ as.list(c(x,params)),
+ rnorm(
+ n=1,
+ mean=n+r*n*(1-n/K)*delta.t,
+ sd=sigma*n*sqrt(delta.t)
+ )
+ )
+ },
+ dprocess=onestep.density,
+ dens.fun=function(x1,x2,t1,t2,params,log,...){
+ delta.t <- t2-t1
+ with(
+ as.list(c(x1,params)),
+ dnorm(
+ x=x2['n'],
+ mean=n+r*n*(1-n/K)*delta.t,
+ sd=sigma*n*sqrt(delta.t),
+ log=log
+ )
+ )
+ },
+ measurement.model=obs~lnorm(meanlog=log(n),sdlog=log(1+tau)),
+ skeleton.vectorfield=function(x,t,params,...){
+ with(
+ as.list(c(x,params)),
+ r*n*(1-n/K)
+ )
+ },
+ delta.t=0.01
+ )
+ simulate(
+ po,
+ params=c(
+ n.0=10000,
+ K=10000,
+ r=0.9,
+ sigma=0.4,
+ tau=0.1
+ ),
+ seed=73658676L
+ )
+ }
+ ) -> verhulst
More information about the pomp-commits
mailing list