[Pomp-commits] r1125 - in pkg/pomp: . src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Mar 2 00:55:22 CET 2015
Author: kingaa
Date: 2015-03-02 00:55:21 +0100 (Mon, 02 Mar 2015)
New Revision: 1125
Added:
pkg/pomp/tests/getting_started.R
pkg/pomp/tests/getting_started.Rout.save
Modified:
pkg/pomp/DESCRIPTION
pkg/pomp/src/pomp_fun.c
Log:
- fix bug causing crashes in rstudio
- add test of 'getting started' vignette's codes
Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION 2015-03-01 16:53:04 UTC (rev 1124)
+++ pkg/pomp/DESCRIPTION 2015-03-01 23:55:21 UTC (rev 1125)
@@ -1,7 +1,7 @@
Package: pomp
Type: Package
Title: Statistical Inference for Partially Observed Markov Processes
-Version: 0.62-3
+Version: 0.62-4
Date: 2015-03-01
Authors at R: c(person(given=c("Aaron","A."),family="King",
role=c("aut","cre"),email="kingaa at umich.edu"),
Modified: pkg/pomp/src/pomp_fun.c
===================================================================
--- pkg/pomp/src/pomp_fun.c 2015-03-01 16:53:04 UTC (rev 1124)
+++ pkg/pomp/src/pomp_fun.c 2015-03-01 23:55:21 UTC (rev 1125)
@@ -36,7 +36,7 @@
}
switch (*mode) {
- case native:
+ case native:
{
SEXP nsi;
PROTECT(nsi = eval(lang3(install("getNativeSymbolInfo"),nf,pack),R_BaseEnv)); nprotect++;
@@ -51,7 +51,7 @@
fname = (const char *) CHARACTER_DATA(STRING_ELT(nf,0));
pkg = (const char *) CHARACTER_DATA(STRING_ELT(pack,0));
ff = R_GetCCallable(pkg,fname);
- PROTECT(f = R_MakeExternalPtr(ff,pfun,R_NilValue)); nprotect++;
+ PROTECT(f = R_MakeExternalPtr(ff,R_NilValue,R_NilValue)); nprotect++;
}
break;
}
Added: pkg/pomp/tests/getting_started.R
===================================================================
--- pkg/pomp/tests/getting_started.R (rev 0)
+++ pkg/pomp/tests/getting_started.R 2015-03-01 23:55:21 UTC (rev 1125)
@@ -0,0 +1,125 @@
+require(pomp)
+
+options(
+ keep.source=TRUE,
+ stringsAsFactors=FALSE,
+ encoding="UTF-8"
+ )
+
+parus.dat <- read.csv(text="
+ year,pop
+ 1960,148
+ 1961,258
+ 1962,185
+ 1963,170
+ 1964,267
+ 1965,239
+ 1966,196
+ 1967,132
+ 1968,167
+ 1969,186
+ 1970,128
+ 1971,227
+ 1972,174
+ 1973,177
+ 1974,137
+ 1975,172
+ 1976,119
+ 1977,226
+ 1978,166
+ 1979,161
+ 1980,199
+ 1981,306
+ 1982,206
+ 1983,350
+ 1984,214
+ 1985,175
+ 1986,211"
+ )
+
+## ----logistic-step-fun---------------------------------------------------
+step.fun <- Csnippet("
+ double dW = rnorm(0,sqrt(dt));
+ N += r*N*(1-N/K)*dt+sigma*N*dW;
+")
+
+## ----logistic-pomp1,cache=FALSE------------------------------------------
+parus <- pomp(data=parus.dat,time="year",t0=1959,
+ rprocess=euler.sim(step.fun=step.fun,delta.t=1/365),
+ statenames="N",paramnames=c("r","K","sigma"))
+
+## ----logistic-simul1-----------------------------------------------------
+simStates <- simulate(parus,nsim=10,params=c(r=0.2,K=200,sigma=0.5,N.0=200),states=TRUE)
+
+## ----logistic-rmeasure---------------------------------------------------
+rmeas <- Csnippet("
+ pop = rpois(phi*N);
+")
+
+## ----logistic-pomp2,cache=FALSE------------------------------------------
+parus <- pomp(parus,rmeasure=rmeas,paramnames="phi",statenames="N")
+
+## ----logistic-simul2-----------------------------------------------------
+sim <- simulate(parus,params=c(r=0.2,K=200,phi=0.5,sigma=0.5,N.0=200),
+ nsim=10,obs=TRUE,states=TRUE)
+
+## ----logistic-dmeasure---------------------------------------------------
+dmeas <- Csnippet("
+ lik = dpois(pop,phi*N,give_log);
+")
+
+## ----logistic-pomp3,cache=FALSE------------------------------------------
+parus <- pomp(parus,dmeasure=dmeas,paramnames="phi",statenames="N")
+
+## ----logistic-pfilter----------------------------------------------------
+pf <- pfilter(parus,Np=1000,params=c(r=0.2,K=200,phi=0.5,sigma=0.5,N.0=200))
+
+## ----logistic-skeleton,cache=FALSE---------------------------------------
+skel <- Csnippet("
+ DN = r*N*(1-N/K);
+")
+
+parus <- pomp(parus,skeleton=skel,skeleton.type="vectorfield",statenames="N",paramnames=c("r","K"))
+
+## ----logistic-traj1------------------------------------------------------
+pars <- parmat(c(r=1,K=200,phi=0.5,sigma=0.5,N.0=20),5)
+pars["N.0",] <- seq(20,300,length=5)
+traj <- trajectory(parus,params=pars,times=seq(1959,1970,by=0.01))
+
+## ----bh-stepfun----------------------------------------------------------
+bh.step <- Csnippet("
+ double eps = rlnorm(-sigma*sigma/2,sigma);
+ N = a*N/(1+b*N)*eps;
+")
+
+## ----bh-skeleton---------------------------------------------------------
+bh.skel <- Csnippet("
+ DN = a*N/(1+b*N);
+")
+
+## ----bh-pomp1,cache=FALSE------------------------------------------------
+parus.bh <- pomp(parus,rprocess=discrete.time.sim(bh.step,delta.t=1),skeleton=bh.skel,skeleton.type="map",skelmap.delta.t=1,statenames="N",paramnames=c("a","b","sigma"))
+
+## ----bh-test-------------------------------------------------------------
+coef(parus.bh) <- c(a=1.1,b=5e-4,sigma=0.5,N.0=30,phi=1)
+sim <- simulate(parus.bh)
+traj <- trajectory(parus.bh)
+pf <- pfilter(parus.bh,Np=1000)
+
+## ----logistic-partrans,cache=FALSE---------------------------------------
+partrans <- Csnippet("
+ Tr = exp(r);
+ TK = exp(K);
+")
+
+parinvtrans <- Csnippet("
+ Tr = log(r);
+ TK = log(K);
+")
+
+parus <- pomp(parus,parameter.transform=partrans,parameter.inv.transform=parinvtrans,paramnames=c("r","K"))
+
+## ----logistic-partrans-test,include=FALSE,cache=FALSE--------------------
+p <- c(r=1,K=200,phi=1,N.0=200,sigma=0.5)
+coef(parus,transform=TRUE) <- partrans(parus,p,dir="inv")
+stopifnot(all.equal(p,coef(parus)))
Added: pkg/pomp/tests/getting_started.Rout.save
===================================================================
--- pkg/pomp/tests/getting_started.Rout.save (rev 0)
+++ pkg/pomp/tests/getting_started.Rout.save 2015-03-01 23:55:21 UTC (rev 1125)
@@ -0,0 +1,149 @@
+
+R version 3.1.2 (2014-10-31) -- "Pumpkin Helmet"
+Copyright (C) 2014 The R Foundation for Statistical Computing
+Platform: x86_64-unknown-linux-gnu (64-bit)
+
+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.
+
+> require(pomp)
+Loading required package: pomp
+Loading required package: subplex
+Loading required package: nloptr
+>
+> options(
++ keep.source=TRUE,
++ stringsAsFactors=FALSE,
++ encoding="UTF-8"
++ )
+>
+> parus.dat <- read.csv(text="
++ year,pop
++ 1960,148
++ 1961,258
++ 1962,185
++ 1963,170
++ 1964,267
++ 1965,239
++ 1966,196
++ 1967,132
++ 1968,167
++ 1969,186
++ 1970,128
++ 1971,227
++ 1972,174
++ 1973,177
++ 1974,137
++ 1975,172
++ 1976,119
++ 1977,226
++ 1978,166
++ 1979,161
++ 1980,199
++ 1981,306
++ 1982,206
++ 1983,350
++ 1984,214
++ 1985,175
++ 1986,211"
++ )
+>
+> ## ----logistic-step-fun---------------------------------------------------
+> step.fun <- Csnippet("
++ double dW = rnorm(0,sqrt(dt));
++ N += r*N*(1-N/K)*dt+sigma*N*dW;
++ ")
+>
+> ## ----logistic-pomp1,cache=FALSE------------------------------------------
+> parus <- pomp(data=parus.dat,time="year",t0=1959,
++ rprocess=euler.sim(step.fun=step.fun,delta.t=1/365),
++ statenames="N",paramnames=c("r","K","sigma"))
+>
+> ## ----logistic-simul1-----------------------------------------------------
+> simStates <- simulate(parus,nsim=10,params=c(r=0.2,K=200,sigma=0.5,N.0=200),states=TRUE)
+>
+> ## ----logistic-rmeasure---------------------------------------------------
+> rmeas <- Csnippet("
++ pop = rpois(phi*N);
++ ")
+>
+> ## ----logistic-pomp2,cache=FALSE------------------------------------------
+> parus <- pomp(parus,rmeasure=rmeas,paramnames="phi",statenames="N")
+>
+> ## ----logistic-simul2-----------------------------------------------------
+> sim <- simulate(parus,params=c(r=0.2,K=200,phi=0.5,sigma=0.5,N.0=200),
++ nsim=10,obs=TRUE,states=TRUE)
+>
+> ## ----logistic-dmeasure---------------------------------------------------
+> dmeas <- Csnippet("
++ lik = dpois(pop,phi*N,give_log);
++ ")
+>
+> ## ----logistic-pomp3,cache=FALSE------------------------------------------
+> parus <- pomp(parus,dmeasure=dmeas,paramnames="phi",statenames="N")
+>
+> ## ----logistic-pfilter----------------------------------------------------
+> pf <- pfilter(parus,Np=1000,params=c(r=0.2,K=200,phi=0.5,sigma=0.5,N.0=200))
+>
+> ## ----logistic-skeleton,cache=FALSE---------------------------------------
+> skel <- Csnippet("
++ DN = r*N*(1-N/K);
++ ")
+>
+> parus <- pomp(parus,skeleton=skel,skeleton.type="vectorfield",statenames="N",paramnames=c("r","K"))
+>
+> ## ----logistic-traj1------------------------------------------------------
+> pars <- parmat(c(r=1,K=200,phi=0.5,sigma=0.5,N.0=20),5)
+> pars["N.0",] <- seq(20,300,length=5)
+> traj <- trajectory(parus,params=pars,times=seq(1959,1970,by=0.01))
+>
+> ## ----bh-stepfun----------------------------------------------------------
+> bh.step <- Csnippet("
++ double eps = rlnorm(-sigma*sigma/2,sigma);
++ N = a*N/(1+b*N)*eps;
++ ")
+>
+> ## ----bh-skeleton---------------------------------------------------------
+> bh.skel <- Csnippet("
++ DN = a*N/(1+b*N);
++ ")
+>
+> ## ----bh-pomp1,cache=FALSE------------------------------------------------
+> parus.bh <- pomp(parus,rprocess=discrete.time.sim(bh.step,delta.t=1),skeleton=bh.skel,skeleton.type="map",skelmap.delta.t=1,statenames="N",paramnames=c("a","b","sigma"))
+>
+> ## ----bh-test-------------------------------------------------------------
+> coef(parus.bh) <- c(a=1.1,b=5e-4,sigma=0.5,N.0=30,phi=1)
+> sim <- simulate(parus.bh)
+> traj <- trajectory(parus.bh)
+> pf <- pfilter(parus.bh,Np=1000)
+>
+> ## ----logistic-partrans,cache=FALSE---------------------------------------
+> partrans <- Csnippet("
++ Tr = exp(r);
++ TK = exp(K);
++ ")
+>
+> parinvtrans <- Csnippet("
++ Tr = log(r);
++ TK = log(K);
++ ")
+>
+> parus <- pomp(parus,parameter.transform=partrans,parameter.inv.transform=parinvtrans,paramnames=c("r","K"))
+>
+> ## ----logistic-partrans-test,include=FALSE,cache=FALSE--------------------
+> p <- c(r=1,K=200,phi=1,N.0=200,sigma=0.5)
+> coef(parus,transform=TRUE) <- partrans(parus,p,dir="inv")
+> stopifnot(all.equal(p,coef(parus)))
+>
+> proc.time()
+ user system elapsed
+ 2.321 0.465 2.752
More information about the pomp-commits
mailing list