[Pomp-commits] r488 - in pkg: data inst/data-R man src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 19 23:50:43 CEST 2011
Author: kingaa
Date: 2011-05-19 23:50:42 +0200 (Thu, 19 May 2011)
New Revision: 488
Modified:
pkg/data/ou2.rda
pkg/inst/data-R/ou2.R
pkg/man/ou2.Rd
pkg/src/ou2.c
pkg/tests/ou2-probe.Rout.save
Log:
- rework the OU2 model slightly
Modified: pkg/data/ou2.rda
===================================================================
(Binary files differ)
Modified: pkg/inst/data-R/ou2.R
===================================================================
--- pkg/inst/data-R/ou2.R 2011-05-19 19:40:43 UTC (rev 487)
+++ pkg/inst/data-R/ou2.R 2011-05-19 21:50:42 UTC (rev 488)
@@ -8,66 +8,13 @@
y2=rep(0,101)
),
t0=1,
- 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 = "ou2_normal_dmeasure",
- rmeasure = "ou2_normal_rmeasure",
+ rprocess=discrete.time.sim("ou2_step",PACKAGE="pomp"),
+ dprocess=onestep.dens("ou2_pdf"),
+ dmeasure = "ou2_dmeasure",
+ rmeasure = "ou2_rmeasure",
skeleton.type="map",
- skeleton = function (x, t, params, ...) {
- with(
- as.list(c(x,params)),
- {
- c(
- x1=alpha.1*x1+alpha.3*x2,
- x2=alpha.2*x1+alpha.4*x2
- )
- }
- )
- },
+ skeleton = "ou2_skel",
+ PACKAGE="pomp",
paramnames = c(
"alpha.1","alpha.2","alpha.3","alpha.4",
"sigma.1","sigma.2","sigma.3",
Modified: pkg/man/ou2.Rd
===================================================================
--- pkg/man/ou2.Rd 2011-05-19 19:40:43 UTC (rev 487)
+++ pkg/man/ou2.Rd 2011-05-19 21:50:42 UTC (rev 488)
@@ -1,30 +1,28 @@
\name{ou2}
\alias{ou2}
\docType{data}
-\title{Two-dimensional Ornstein-Uhlenbeck process}
+\title{Two-dimensional discrete-time Ornstein-Uhlenbeck process}
\description{
\code{ou2} is a \code{pomp} object encoding a bivariate discrete-time Ornstein-Uhlenbeck process.
}
\usage{data(ou2)}
\details{
- The process is fully but noisily observed.
+ If the state process is \eqn{X(t) = (x_{1}(t),x_{2}(t))}, then
+ \deqn{X(t+1) = \alpha X(t) + \sigma \epsilon(t),}
+ where \eqn{\alpha} and \eqn{\sigma} are 2x2 matrices, \eqn{\sigma} is lower-triangular, and \eqn{\epsilon(t)} is standard bivariate normal.
+ The observation process is \eqn{Y(t) = (y_1(t),y_2(t))}, where \eqn{y_i(t) \sim \mathrm{normal}(x_i(t),\tau)}.
The functions \code{rprocess}, \code{dprocess}, \code{rmeasure}, \code{dmeasure}, and \code{skeleton} are implemented using compiled C code for computational speed:
see the source code for details.
- The use of this object is demonstrated in the vignette "intro_to_pomp".
+ This object is demonstrated in the vignette "Advanced topics in pomp".
}
\examples{
data(ou2)
plot(ou2)
coef(ou2)
-theta <- c(
- alpha.1=0.8, alpha.2=-0.5, alpha.3=0.3, alpha.4=0.9,
- sigma.1=3, sigma.2=-0.5, sigma.3=2,
- tau=1,
- x1.0=-3, x2.0=4
- )
-x <- simulate(ou2,nsim=10,seed=20348585,params=theta)
-plot(x[[1]])
-pfilter(ou2,params=theta,Np=1000)$loglik
+x <- simulate(ou2)
+plot(x)
+pf <- pfilter(ou2,Np=1000)
+logLik(pf)
}
\seealso{\code{\link{pomp}} and the vignettes}
\keyword{datasets}
Modified: pkg/src/ou2.c
===================================================================
--- pkg/src/ou2.c 2011-05-19 19:40:43 UTC (rev 487)
+++ pkg/src/ou2.c 2011-05-19 21:50:42 UTC (rev 488)
@@ -5,23 +5,25 @@
#include <Rdefines.h>
// prototypes
-void ou2_normal_rmeasure (double *y, double *x, double *p,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
- int ncovar, double *covar, double t);
-void ou2_normal_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
- int covdim, double *covar, double t);
+void ou2_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovar, double *covar, double t);
+void ou2_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int covdim, double *covar, double t);
void ou2_adv (double *x, double *xstart, double *par, double *times, int *n);
-void _ou2_adv (double *x, double *xstart, double *par, double *times, int *n, int *parindex);
void ou2_step (double *x, const double *p,
const int *stateindex, const int *parindex, const int *covindex,
int ncovars, const double *covars,
double t, double dt);
-void ou2_pdf (double *d, double *X, double *par, double *times, int *n, int *parindex, int *give_log);
-static void sim_ou2 (double *x,
+void ou2_pdf (double *f,
+ double *x1, double *x2, double t1, double t2, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int ncovars, const double *covars);
+static void sim_ou2 (double *x1, double *x2,
double alpha1, double alpha2, double alpha3, double alpha4,
double sigma1, double sigma2, double sigma3);
-static double dens_ou2 (double *x1, double *x2,
+static double dens_ou2 (double x1, double x2, double z1, double z2,
double alpha1, double alpha2, double alpha3, double alpha4,
double sigma1, double sigma2, double sigma3, int give_log);
@@ -70,7 +72,7 @@
tgoal = times[k];
while (tnow < tgoal) {
- sim_ou2(xp1,ALPHA1,ALPHA2,ALPHA3,ALPHA4,SIGMA1,SIGMA2,SIGMA3); // advance state
+ sim_ou2(&xp1[0],&xp1[1],ALPHA1,ALPHA2,ALPHA3,ALPHA4,SIGMA1,SIGMA2,SIGMA3); // advance state
tnow += dt; // advance time
}
@@ -89,85 +91,6 @@
#undef SIGMA2
#undef SIGMA3
-#define ALPHA1 (pp[parindex[0]])
-#define ALPHA2 (pp[parindex[1]])
-#define ALPHA3 (pp[parindex[2]])
-#define ALPHA4 (pp[parindex[3]])
-#define SIGMA1 (pp[parindex[4]])
-#define SIGMA2 (pp[parindex[5]])
-#define SIGMA3 (pp[parindex[6]])
-
-// just like the above, but with parindex
-void _ou2_adv (double *x, double *xstart, double *par, double *times, int *n, int *parindex)
-{
- int nvar = n[0], npar = n[1], nrep = n[2], ntimes = n[3], incr;
- double tnow, tgoal, dt = 1.0;
- double *xp0, *xp1, *pp;
- int i, j, k;
-
- incr = nrep*nvar;
-
- GetRNGstate(); // initialize R's pseudorandom number generator
-
- for (j = 0; j < nrep; j++) {
-
- R_CheckUserInterrupt(); // check for an interrupt signal
-
- xp0 = &xstart[nvar*j]; // pointer to j-th starting state
- xp1 = &x[nvar*j]; // pointer to j-th state vector
- pp = &par[npar*j]; // pointer to j-th parameter vector
-
- for (i = 0; i < nvar; i++) xp1[i] = xp0[i]; // copy xstart into the first slice of x
-
- tnow = times[0]; // initial time
-
- for (k = 1; k < ntimes; k++) { // loop over times
-
- xp0 = xp1;
- xp1 += incr;
-
- for (i = 0; i < nvar; i++) xp1[i] = xp0[i]; // copy state vector
-
- tgoal = times[k];
-
- while (tnow < tgoal) {
- sim_ou2(xp1,ALPHA1,ALPHA2,ALPHA3,ALPHA4,SIGMA1,SIGMA2,SIGMA3); // advance state
- tnow += dt; // advance time
- }
-
- }
-
- }
-
- PutRNGstate(); // finished with R's random number generator
-}
-
-// pdf of a single 2D OU transition
-void ou2_pdf (double *d, double *X, double *par, double *times, int *n, int *parindex, int *give_log)
-{
- int nvar = n[0], npar = n[1], nrep = n[2], ntimes = n[3];
- double *x1, *x2, *pp;
- int j, k;
- for (k = 0; k < nrep; k++) {
- pp = &par[npar*k]; // get address of k-th parameter vector
- x1 = &X[nvar*k]; // get address of (0,0)-th state vector
- for (j = 1; j < ntimes; j++) {
- R_CheckUserInterrupt();
- x2 = &X[nvar*(k+nrep*j)]; // get address of (k,j)-th state vector
- d[k+nrep*(j-1)] = dens_ou2(x1,x2,ALPHA1,ALPHA2,ALPHA3,ALPHA4,SIGMA1,SIGMA2,SIGMA3,*give_log);
- x1 = x2;
- }
- }
-}
-
-#undef ALPHA1
-#undef ALPHA2
-#undef ALPHA3
-#undef ALPHA4
-#undef SIGMA1
-#undef SIGMA2
-#undef SIGMA3
-
#define ALPHA1 (p[parindex[0]])
#define ALPHA2 (p[parindex[1]])
#define ALPHA3 (p[parindex[2]])
@@ -177,8 +100,9 @@
#define SIGMA3 (p[parindex[6]])
#define TAU (p[parindex[7]])
-#define X1 (x[stateindex[0]])
-#define X2 (x[stateindex[1]])
+#define X1 (stateindex[0])
+#define X2 (stateindex[1])
+
#define Y1 (y[obsindex[0]])
#define Y2 (y[obsindex[1]])
@@ -188,31 +112,50 @@
int ncovars, const double *covars,
double t, double dt)
{
- sim_ou2(x,ALPHA1,ALPHA2,ALPHA3,ALPHA4,SIGMA1,SIGMA2,SIGMA3);
+ sim_ou2(&x[X1],&x[X2],ALPHA1,ALPHA2,ALPHA3,ALPHA4,SIGMA1,SIGMA2,SIGMA3);
}
+// onestep transition probability density for use in 'onestep.dens' plug-in
+// transition from x to z as time goes from t1 to t2
+void ou2_pdf (double *f,
+ double *x, double *z, double t1, double t2, const double *p,
+ const int *stateindex, const int *parindex, const int *covindex,
+ int ncovars, const double *covars)
+{
+ if (t2-t1 != 1)
+ error("ou2_pdf error: transitions must be consecutive");
+ f[0] = dens_ou2(x[X1],x[X2],z[X1],z[X2],ALPHA1,ALPHA2,ALPHA3,ALPHA4,SIGMA1,SIGMA2,SIGMA3,1);
+}
+void ou2_skel (double *f, double *x, double *p,
+ int *stateindex, int *parindex, int *covindex,
+ int ncovars, double *covars, double t)
+{
+ f[X1] = ALPHA1*x[X1]+ALPHA3*x[X2];
+ f[X2] = ALPHA2*x[X1]+ALPHA4*x[X2];
+}
+
// bivariate normal measurement error density
-void ou2_normal_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
- int covdim, double *covar, double t)
+void ou2_dmeasure (double *lik, double *y, double *x, double *p, int give_log,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int covdim, double *covar, double t)
{
double sd = fabs(TAU);
double f = 0.0;
- f += (ISNA(Y1)) ? 0.0 : dnorm(Y1,X1,sd,1);
- f += (ISNA(Y2)) ? 0.0 : dnorm(Y2,X2,sd,1);
+ f += (ISNA(Y1)) ? 0.0 : dnorm(Y1,x[X1],sd,1);
+ f += (ISNA(Y2)) ? 0.0 : dnorm(Y2,x[X2],sd,1);
*lik = (give_log) ? f : exp(f);
}
// bivariate normal measurement error simulator
-void ou2_normal_rmeasure (double *y, double *x, double *p,
- int *obsindex, int *stateindex, int *parindex, int *covindex,
- int ncovar, double *covar,
- double t)
+void ou2_rmeasure (double *y, double *x, double *p,
+ int *obsindex, int *stateindex, int *parindex, int *covindex,
+ int ncovar, double *covar,
+ double t)
{
double sd = fabs(TAU);
- Y1 = rnorm(X1,sd);
- Y2 = rnorm(X2,sd);
+ Y1 = rnorm(x[X1],sd);
+ Y2 = rnorm(x[X2],sd);
}
#undef ALPHA1
@@ -230,14 +173,14 @@
#undef Y2
// simple 2D Ornstein-Uhlenbeck process simulation
-static void sim_ou2 (double *x,
+static void sim_ou2 (double *x1, double *x2,
double alpha1, double alpha2, double alpha3, double alpha4,
double sigma1, double sigma2, double sigma3)
{
double eps[2], xnew[2];
- if (!(R_FINITE(x[0]))) return;
- if (!(R_FINITE(x[1]))) return;
+ if (!(R_FINITE(*x1))) return;
+ if (!(R_FINITE(*x2))) return;
if (!(R_FINITE(alpha1))) return;
if (!(R_FINITE(alpha2))) return;
if (!(R_FINITE(alpha3))) return;
@@ -249,24 +192,25 @@
eps[0] = rnorm(0,1);
eps[1] = rnorm(0,1);
- xnew[0] = alpha1*x[0]+alpha3*x[1]+sigma1*eps[0];
- xnew[1] = alpha2*x[0]+alpha4*x[1]+sigma2*eps[0]+sigma3*eps[1];
+ xnew[0] = alpha1*(*x1)+alpha3*(*x2)+sigma1*eps[0];
+ xnew[1] = alpha2*(*x1)+alpha4*(*x2)+sigma2*eps[0]+sigma3*eps[1];
- x[0] = xnew[0];
- x[1] = xnew[1];
+ *x1 = xnew[0];
+ *x2 = xnew[1];
}
-// simple 2D Ornstein-Uhlenbeck process density
-static double dens_ou2 (double *x1, double *x2,
+// simple 2D Ornstein-Uhlenbeck process transition density
+// transition (x1,x2) -> (z1,z2) in 1 unit of time
+static double dens_ou2 (double x1, double x2, double z1, double z2,
double alpha1, double alpha2, double alpha3, double alpha4,
double sigma1, double sigma2, double sigma3, int give_log)
{
double eps[2], val;
- if (!(R_FINITE(x1[0]))) return R_NaReal;
- if (!(R_FINITE(x1[1]))) return R_NaReal;
- if (!(R_FINITE(x2[0]))) return R_NaReal;
- if (!(R_FINITE(x2[1]))) return R_NaReal;
+ if (!(R_FINITE(x1))) return R_NaReal;
+ if (!(R_FINITE(x2))) return R_NaReal;
+ if (!(R_FINITE(z1))) return R_NaReal;
+ if (!(R_FINITE(z2))) return R_NaReal;
if (!(R_FINITE(alpha1))) return R_NaReal;
if (!(R_FINITE(alpha2))) return R_NaReal;
if (!(R_FINITE(alpha3))) return R_NaReal;
@@ -276,8 +220,8 @@
if (!(R_FINITE(sigma3))) return R_NaReal;
// compute residuals
- eps[0] = x2[0]-alpha1*x1[0]-alpha3*x1[1];
- eps[1] = x2[1]-alpha2*x1[0]-alpha4*x1[1];
+ eps[0] = z1-alpha1*x1-alpha3*x2;
+ eps[1] = z2-alpha2*x1-alpha4*x2;
// backsolve
eps[0] /= sigma1;
Modified: pkg/tests/ou2-probe.Rout.save
===================================================================
--- pkg/tests/ou2-probe.Rout.save 2011-05-19 19:40:43 UTC (rev 487)
+++ pkg/tests/ou2-probe.Rout.save 2011-05-19 21:50:42 UTC (rev 488)
@@ -62,14 +62,14 @@
$quantiles
y1.mean y2.mean y1.sd y2.sd
- 0.122 0.534 0.028 0.042
+ 0.126 0.518 0.034 0.076
$pvals
y1.mean y2.mean y1.sd y2.sd
-0.24750499 0.93413174 0.05988024 0.08782435
+0.25548902 0.96606786 0.07185629 0.15568862
$synth.loglik
-[1] -4.780494
+[1] -4.347702
> summary(pm.po)
$coef
@@ -81,14 +81,14 @@
$quantiles
y1.mean y2.mean y1.sd y2.sd
- 0.078 0.682 1.000 1.000
+ 0.076 0.662 1.000 1.000
$pvals
y1.mean y2.mean y1.sd y2.sd
-0.159680639 0.638722555 0.003992016 0.003992016
+0.155688623 0.678642715 0.003992016 0.003992016
$synth.loglik
-[1] -111.0578
+[1] -108.7680
>
> plot(pm.ou2)
@@ -117,18 +117,18 @@
$quantiles
y1acf.acf.1.y1 y1acf.acf.2.y1 y2acf.acf.0.y2 y2acf.acf.1.y2 y2acf.acf.2.y2
- 0.098 0.040 0.066 0.068 0.066
+ 0.106 0.040 0.062 0.066 0.068
y12ccf.ccf.3 y12ccf.ccf.8
- 0.052 0.324
+ 0.054 0.362
$pvals
y1acf.acf.1.y1 y1acf.acf.2.y1 y2acf.acf.0.y2 y2acf.acf.1.y2 y2acf.acf.2.y2
- 0.19960080 0.08383234 0.13572854 0.13972056 0.13572854
+ 0.21556886 0.08383234 0.12774451 0.13572854 0.13972056
y12ccf.ccf.3 y12ccf.ccf.8
- 0.10778443 0.65069860
+ 0.11177645 0.72654691
$synth.loglik
-[1] -9.00962
+[1] -9.302747
>
> pb <- probe(
@@ -150,22 +150,22 @@
$quantiles
y1.10% y1.20% y1.30% y1.40% y1.50% y1.60% y1.70% y1.80%
- 0.730 0.770 0.770 0.545 0.355 0.135 0.020 0.010
+ 0.735 0.750 0.810 0.585 0.325 0.135 0.025 0.005
y1.90% acf.0.y1 acf.1.y1 acf.4.y1 acf.7.y1 acf.0.y2 acf.1.y2 acf.4.y2
- 0.000 0.030 0.035 0.835 0.860 0.070 0.075 0.465
+ 0.000 0.030 0.050 0.815 0.855 0.075 0.070 0.420
acf.7.y2 pd
- 0.840 0.100
+ 0.870 0.080
$pvals
y1.10% y1.20% y1.30% y1.40% y1.50% y1.60%
-0.547263682 0.467661692 0.467661692 0.915422886 0.716417910 0.278606965
+0.537313433 0.507462687 0.388059701 0.835820896 0.656716418 0.278606965
y1.70% y1.80% y1.90% acf.0.y1 acf.1.y1 acf.4.y1
-0.049751244 0.029850746 0.009950249 0.069651741 0.079601990 0.338308458
+0.059701493 0.019900498 0.009950249 0.069651741 0.109452736 0.378109453
acf.7.y1 acf.0.y2 acf.1.y2 acf.4.y2 acf.7.y2 pd
-0.288557214 0.149253731 0.159203980 0.935323383 0.328358209 0.208955224
+0.298507463 0.159203980 0.149253731 0.845771144 0.268656716 0.169154229
$synth.loglik
-[1] -26.82926
+[1] -25.89253
>
> po <- ou2
@@ -216,14 +216,14 @@
$quantiles
acf.0.y1 acf.1.y1 acf.2.y1 ccf.0 ccf.1 ccf.2
- 0.850 0.853 0.840 0.850 0.853 0.840
+ 0.860 0.857 0.840 0.860 0.857 0.840
$pvals
acf.0.y1 acf.1.y1 acf.2.y1 ccf.0 ccf.1 ccf.2
-0.3016983 0.2957043 0.3216783 0.3016983 0.2957043 0.3216783
+0.2817183 0.2877123 0.3216783 0.2817183 0.2877123 0.3216783
$synth.loglik
-[1] 89.34863
+[1] 91.41242
>
> pb <- probe(
@@ -243,14 +243,14 @@
$quantiles
ccf.-5 ccf.-3 ccf.1 ccf.4 ccf.8
- 0.834 0.860 0.149 0.164 0.929
+ 0.838 0.860 0.123 0.157 0.912
$pvals
ccf.-5 ccf.-3 ccf.1 ccf.4 ccf.8
-0.3336663 0.2817183 0.2997003 0.3296703 0.1438561
+0.3256743 0.2817183 0.2477522 0.3156843 0.1778222
$synth.loglik
-[1] -14.9071
+[1] -14.50221
>
> pb <- probe(
@@ -270,24 +270,24 @@
$quantiles
ccf.-5 ccf.-3 ccf.1 ccf.4 ccf.8
- 0.693 0.821 0.346 0.243 0.850
+ 0.728 0.819 0.343 0.238 0.840
$pvals
ccf.-5 ccf.-3 ccf.1 ccf.4 ccf.8
-0.6153846 0.3596404 0.6933067 0.4875125 0.3016983
+0.5454545 0.3636364 0.6873127 0.4775225 0.3216783
$synth.loglik
-[1] 7.363337
+[1] 7.505796
>
> head(as(pb,"data.frame"))
- ccf.-5 ccf.-3 ccf.1 ccf.4 ccf.8
-data 0.6340061 0.8339776 -0.3925899 -0.7963667 0.301628540
-sim.1 0.4470813 0.8292044 -0.3798152 -0.7420235 0.519633798
-sim.2 0.5461973 0.7227145 -0.2996333 -0.6294983 0.012241523
-sim.3 0.4229191 0.7283615 -0.1026817 -0.5565477 -0.115334617
-sim.4 0.4577827 0.7103028 -0.3372048 -0.5293968 -0.020542485
-sim.5 0.4013560 0.6594786 -0.2933890 -0.5947615 -0.001712792
+ ccf.-5 ccf.-3 ccf.1 ccf.4 ccf.8
+data 0.6340061 0.8339776 -0.3925899 -0.7963667 0.30162854
+sim.1 0.4727712 0.8457682 -0.3980453 -0.7878552 0.42505373
+sim.2 0.6416744 0.8419569 -0.3851185 -0.7798302 0.25492015
+sim.3 0.6419382 0.8434833 -0.3211908 -0.7538791 0.15241925
+sim.4 0.3693204 0.5205152 -0.4305519 -0.5773910 -0.04953301
+sim.5 0.5542739 0.6216561 -0.4089672 -0.7236761 0.06850115
>
>
>
More information about the pomp-commits
mailing list