[Pomp-commits] r1027 - in pkg/pompExamples: . inst/examples src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Dec 17 20:12:41 CET 2014
Author: kingaa
Date: 2014-12-17 20:12:41 +0100 (Wed, 17 Dec 2014)
New Revision: 1027
Added:
pkg/pompExamples/src/bbp.c
Modified:
pkg/pompExamples/DESCRIPTION
pkg/pompExamples/inst/examples/bbp.R
pkg/pompExamples/tests/examples.R
pkg/pompExamples/tests/pertussis.Rout.save
Log:
- work on Bombay plague example
Modified: pkg/pompExamples/DESCRIPTION
===================================================================
--- pkg/pompExamples/DESCRIPTION 2014-12-17 19:12:36 UTC (rev 1026)
+++ pkg/pompExamples/DESCRIPTION 2014-12-17 19:12:41 UTC (rev 1027)
@@ -17,7 +17,6 @@
URL: http://pomp.r-forge.r-project.org
Description: More 'pomp' examples.
Depends: R(>= 3.0.0), stats, graphics, pomp(>= 0.49-1)
-Suggests: Rmpi,mpifarm,plyr,reshape2,ggplot2,xtable
License: GPL (>= 2)
LazyData: false
BuildVignettes: true
Modified: pkg/pompExamples/inst/examples/bbp.R
===================================================================
--- pkg/pompExamples/inst/examples/bbp.R 2014-12-17 19:12:36 UTC (rev 1026)
+++ pkg/pompExamples/inst/examples/bbp.R 2014-12-17 19:12:41 UTC (rev 1027)
@@ -63,12 +63,23 @@
y += dy + beta*X*(dW+beta*X*ito);
n += dn;
"
+ ),
+ delta.t=1/24/7
),
- delta.t=1/24),
+ skeleton=Csnippet("
+ double X = exp(x);
+ double Y = exp(y);
+ Dx = mu*(1.0/X-1)+(delta-beta)*Y;
+ Dy = beta*X+delta*(Y-1)-gamma-mu;
+ Dn = -delta*Y;
+ "
+ ),
+ skeleton.type="vectorfield",
paramnames=c("beta","delta","mu","gamma","sigma","theta","ratio"),
statenames=c("x","y","n"),
- measurement.model=deaths~nbinom(mu=ratio*exp(y),size=theta),
- logvar=c("beta","delta","ratio","sigma","theta"),
+ rmeasure=Csnippet("deaths=rnbinom_mu(theta,ratio*exp(y));"),
+ dmeasure=Csnippet("lik=dnbinom_mu(deaths,theta,ratio*exp(y),give_log);"),
+ logvar=c("beta","delta","ratio","sigma","theta","mu"),
logitvar=c("y0"),
parameter.inv.transform=function (params, logvar, logitvar, ...) {
params[logvar] <- log(params[logvar])
Added: pkg/pompExamples/src/bbp.c
===================================================================
--- pkg/pompExamples/src/bbp.c (rev 0)
+++ pkg/pompExamples/src/bbp.c 2014-12-17 19:12:41 UTC (rev 1027)
@@ -0,0 +1,103 @@
+/* pomp model file: _bombay_plague */
+
+#include <pomp.h>
+#include <R_ext/Rdynload.h>
+
+
+#define Beta (__p[__parindex[0]])
+#define delta (__p[__parindex[1]])
+#define mu (__p[__parindex[2]])
+#define gamma (__p[__parindex[3]])
+#define sigma (__p[__parindex[4]])
+#define theta (__p[__parindex[5]])
+#define ratio (__p[__parindex[6]])
+#define x (__x[__stateindex[0]])
+#define y (__x[__stateindex[1]])
+#define n (__x[__stateindex[2]])
+#define deaths (__y[__obsindex[0]])
+#define Dx (__f[__stateindex[0]])
+#define Dy (__f[__stateindex[1]])
+#define Dn (__f[__stateindex[2]])
+#define TBeta (__pt[__parindex[0]])
+#define Tdelta (__pt[__parindex[1]])
+#define Tmu (__pt[__parindex[2]])
+#define Tgamma (__pt[__parindex[3]])
+#define Tsigma (__pt[__parindex[4]])
+#define Ttheta (__pt[__parindex[5]])
+#define Tratio (__pt[__parindex[6]])
+#define lik (__lik[0])
+
+void _bombay_plague_rmeasure (double *__y, double *__x, double *__p, int *__obsindex, int *__stateindex, int *__parindex, int *__covindex, int __ncovars, double *__covars, double t)
+{
+ deaths=rnbinom_mu(theta,ratio*exp(y));
+}
+
+
+void _bombay_plague_dmeasure (double *__lik, double *__y, double *__x, double *__p, int give_log, int *__obsindex, int *__stateindex, int *__parindex, int *__covindex, int __ncovars, double *__covars, double t)
+{
+ lik=dnbinom_mu(deaths,theta,ratio*exp(y),give_log);
+}
+
+
+void _bombay_plague_stepfn (double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, int __covdim, const double *__covars, double t, double dt)
+{
+
+ double X = exp(x);
+ double Y = exp(y);
+ double dx, dy, dn, dW, ito;
+ dx = (mu*(1.0/X-1)+(delta-Beta)*Y)*dt;
+ dy = (Beta*X+delta*(Y-1)-gamma-mu)*dt;
+ dn = -delta*Y*dt;
+ dW = rnorm(0,sigma*sqrt(dt));
+ ito = 0.5*sigma*sigma*dt;
+ x += dx - Beta*Y*(dW-Beta*Y*ito);
+ y += dy + Beta*X*(dW+Beta*X*ito);
+ n += dn;
+
+}
+
+
+void _bombay_plague_skelfn (double *__f, double *__x, double *__p, int *__stateindex, int *__parindex, int *__covindex, int __ncovars, double *__covars, double t)
+{
+
+ double X = exp(x);
+ double Y = exp(y);
+ Dx = mu*(1.0/X-1)+(delta-Beta)*Y;
+ Dy = Beta*X+delta*(Y-1)-gamma-mu;
+ Dn = -delta*Y;
+
+}
+
+
+void _bombay_plague_rprior (double *__p, int *__parindex)
+{
+ error("'rprior' not defined");
+}
+
+
+void _bombay_plague_dprior (double *__lik, double *__p, int give_log, int *__parindex)
+{
+ error("'dprior' not defined");
+}
+
+#undef Beta
+#undef delta
+#undef mu
+#undef gamma
+#undef sigma
+#undef theta
+#undef ratio
+#undef x
+#undef y
+#undef n
+#undef deaths
+#undef Dx
+#undef Dy
+#undef Dn
+#undef TBeta
+#undef Tdelta
+#undef Tmu
+#undef Tgamma
+#undef Tsigma
+#undef Ttheta
+#undef Tratio
Modified: pkg/pompExamples/tests/examples.R
===================================================================
--- pkg/pompExamples/tests/examples.R 2014-12-17 19:12:36 UTC (rev 1026)
+++ pkg/pompExamples/tests/examples.R 2014-12-17 19:12:41 UTC (rev 1027)
@@ -21,3 +21,4 @@
pompExample(bbp)
pf <- pfilter(simulate(bbp),Np=100,max.fail=Inf)
+tj <- trajectory(bbp)
Modified: pkg/pompExamples/tests/pertussis.Rout.save
===================================================================
--- pkg/pompExamples/tests/pertussis.Rout.save 2014-12-17 19:12:36 UTC (rev 1026)
+++ pkg/pompExamples/tests/pertussis.Rout.save 2014-12-17 19:12:41 UTC (rev 1027)
@@ -1,5 +1,5 @@
-R version 3.1.2 (2014-10-31) -- "Pumpkin Helmet"
+R Under development (unstable) (2014-12-14 r67168) -- "Unsuffered Consequences"
Copyright (C) 2014 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
@@ -146,7 +146,7 @@
>
> system.time(pf <- pfilter(pertussis.sim(full.small),seed=3445886L,Np=1000))
user system elapsed
- 18.793 0.000 18.864
+ 17.801 0.004 17.861
> logLik(pf)
[1] -3829.33
>
@@ -170,4 +170,4 @@
>
> proc.time()
user system elapsed
- 19.533 0.040 19.671
+ 18.657 0.060 18.801
More information about the pomp-commits
mailing list