[Pomp-commits] r389 - in pkg: R man tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Oct 19 15:49:03 CEST 2010
Author: kingaa
Date: 2010-10-19 15:49:03 +0200 (Tue, 19 Oct 2010)
New Revision: 389
Modified:
pkg/R/probe.R
pkg/man/probed-pomp-class.Rd
pkg/tests/ou2-probe.Rout.save
pkg/tests/ricker-probe.Rout.save
Log:
- add 'synth.loglik' slot to the 'probed.pomp' object
Modified: pkg/R/probe.R
===================================================================
--- pkg/R/probe.R 2010-10-19 13:21:24 UTC (rev 388)
+++ pkg/R/probe.R 2010-10-19 13:49:03 UTC (rev 389)
@@ -7,7 +7,8 @@
datvals="numeric",
simvals="array",
quantiles="numeric",
- pvals="numeric"
+ pvals="numeric",
+ synth.loglik="numeric"
)
)
@@ -19,7 +20,8 @@
coef=coef(object),
nsim=nrow(object at simvals),
quantiles=object at quantiles,
- pvals=object at pvals
+ pvals=object at pvals,
+ synth.loglik=object at synth.loglik
)
}
)
@@ -70,6 +72,8 @@
quants[k] <- sum(simval[,k]<datval[k])/nsim
}
+ ll <- .Call(synth_loglik,simval,datval)
+
coef(object) <- params
new(
@@ -80,7 +84,8 @@
datvals=datval,
simvals=simval,
quantiles=quants,
- pvals=pvals
+ pvals=pvals,
+ synth.loglik=ll
)
}
)
Modified: pkg/man/probed-pomp-class.Rd
===================================================================
--- pkg/man/probed-pomp-class.Rd 2010-10-19 13:21:24 UTC (rev 388)
+++ pkg/man/probed-pomp-class.Rd 2010-10-19 13:49:03 UTC (rev 389)
@@ -22,6 +22,7 @@
\item{datvals, simvals}{values of each of the probes applied to the real and simulated data, respectively.}
\item{quantiles}{fraction of simulations with probe values less than the value of the probe of the data.}
\item{pvals}{two-sided p-values: fraction of the \code{simvals} that deviate more extremely from the mean of the \code{simvals} than does \code{datavals}.}
+ \item{synth.loglik}{the log synthetic likelihood (Wood 2010). This is the likelihood assuming that the probes are multivariate-normally distributed.}
\item{weights}{relative weights applied to the discrepancies in computing the probe-matching objective function.}
\item{fail.value}{value to use to replace non-finite values of the objective function.}
\item{evals}{
@@ -52,6 +53,10 @@
B. E. Kendall, C. J. Briggs, W. M. Murdoch, P. Turchin, S. P. Ellner, E. McCauley, R. M. Nisbet, S. N. Wood
Why do populations cycle? A synthesis of statistical and mechanistic modeling approaches,
Ecology, 80:1789--1805, 1999.
+
+ S. N. Wood
+ Statistical inference for noisy nonlinear ecological dynamic systems,
+ Nature, 466: 1102--1104, 2010.
}
\author{
Daniel C. Reuman \email{d.reuman at imperial dot ac dot uk}
Modified: pkg/tests/ou2-probe.Rout.save
===================================================================
--- pkg/tests/ou2-probe.Rout.save 2010-10-19 13:21:24 UTC (rev 388)
+++ pkg/tests/ou2-probe.Rout.save 2010-10-19 13:49:03 UTC (rev 389)
@@ -30,6 +30,7 @@
+ ),
+ nsim=500
+ )
+>
> pm.po <- probe(
+ ou2,
+ params=c(
@@ -63,6 +64,9 @@
y1.mean y2.mean y1.sd y2.sd
0.24750499 0.91816367 0.04391218 0.10379242
+$synth.loglik
+[1] -5.480174
+
> summary(pm.po)
$coef
alpha.1 alpha.2 alpha.3 alpha.4 sigma.1 sigma.2 sigma.3 tau x1.0 x2.0
@@ -79,6 +83,9 @@
y1.mean y2.mean y1.sd y2.sd
0.151696607 0.702594810 0.003992016 0.003992016
+$synth.loglik
+[1] -180.4989
+
>
> plot(pm.ou2)
> plot(pm.po)
@@ -116,6 +123,9 @@
y12ccf.ccf.3 y12ccf.ccf.8
0.09980040 0.65469062
+$synth.loglik
+[1] -3368.175
+
>
> pb <- probe(
+ ou2,
@@ -150,6 +160,9 @@
acf.1.y2 acf.4.y2 acf.7.y2 pd
0.14925373 0.86567164 0.26865672 0.12935323
+$synth.loglik
+[1] -1.732555e+39
+
>
> po <- ou2
> coef(po,c("alpha.2","alpha.3")) <- c(0,0)
@@ -183,7 +196,7 @@
+ po,
+ probes=list(
+ probe.acf(var=c("y1"),lags=c(0,1,2),type="cov"),
-+ probe.ccf(vars=c("y1","y1"),lags=c(0,1,2))
++ probe.ccf(vars=c("y1","y1"),lags=c(0,1,2),type="cov")
+ ),
+ nsim=1000,
+ seed=1066L
@@ -205,6 +218,9 @@
acf.0.y1 acf.1.y1 acf.2.y1 ccf.0 ccf.1 ccf.2
0.1918082 0.2237762 0.4495504 0.1918082 0.2237762 0.4495504
+$synth.loglik
+[1] 27.6966
+
>
> pb <- probe(
+ po,
@@ -229,15 +245,45 @@
ccf.-5 ccf.-3 ccf.1 ccf.4 ccf.8
0.4075924 0.3876124 0.0919081 0.3856144 0.3716284
+$synth.loglik
+[1] -27.84814
+
>
+> pb <- probe(
++ po,
++ probes=probe.ccf(vars=c("y1","y2"),lags=c(-5,-3,1,4,8),type="corr"),
++ nsim=1000,
++ seed=1066L
++ )
+> plot(pb)
+> summary(pb)
+$coef
+alpha.1 alpha.2 alpha.3 alpha.4 sigma.1 sigma.2 sigma.3 tau x1.0 x2.0
+ 0.8 -0.5 0.3 0.9 3.0 -0.5 2.0 1.0 -3.0 4.0
+
+$nsim
+[1] 1000
+
+$quantiles
+ccf.-5 ccf.-3 ccf.1 ccf.4 ccf.8
+ 0.364 0.425 0.961 0.671 0.176
+
+$pvals
+ ccf.-5 ccf.-3 ccf.1 ccf.4 ccf.8
+0.72927073 0.85114885 0.07992008 0.65934066 0.35364635
+
+$synth.loglik
+[1] 2.860594
+
+>
> head(as(pb,"data.frame"))
- ccf.-5 ccf.-3 ccf.1 ccf.4 ccf.8
-data 19.34604 27.20885 -9.222381 -24.94993 -0.74878977
-sim.1 25.39116 47.09313 -21.570904 -42.14185 29.51164251
-sim.2 20.36428 26.94550 -11.171452 -23.47005 0.45640982
-sim.3 11.55157 19.89440 -2.804638 -15.20150 -3.15023990
-sim.4 21.26787 32.99956 -15.666009 -24.59495 -0.95437169
-sim.5 13.42573 22.06017 -9.814133 -19.89532 -0.05729449
+ ccf.-5 ccf.-3 ccf.1 ccf.4 ccf.8
+data 0.5333071 0.7500589 -0.2542309 -0.6877880 -0.020641683
+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
>
>
>
Modified: pkg/tests/ricker-probe.Rout.save
===================================================================
--- pkg/tests/ricker-probe.Rout.save 2010-10-19 13:21:24 UTC (rev 388)
+++ pkg/tests/ricker-probe.Rout.save 2010-10-19 13:49:03 UTC (rev 389)
@@ -47,6 +47,9 @@
$pvals
[1] 0.3496503
+$synth.loglik
+[1] -2.437844
+
>
> pb <- probe(
+ po,
@@ -71,6 +74,9 @@
nlar.1^1 nlar.2^1 nlar.3^1
0.3336663 0.5034965 0.3956044
+$synth.loglik
+[1] 0.2080522
+
>
> pb <- probe(
+ po,
@@ -95,6 +101,9 @@
nlar.1^1 nlar.2^1 nlar.3^1
0.3336663 0.5034965 0.3956044
+$synth.loglik
+[1] 0.2080522
+
>
> pb <- probe(
+ po,
@@ -119,6 +128,9 @@
nlar.1^1 nlar.1^2 nlar.1^3
0.8071928 0.3916084 0.6473526
+$synth.loglik
+[1] 9.406785
+
>
> pb <- probe(
+ po,
@@ -151,6 +163,9 @@
marg.1 marg.2 marg.3
0.2917083 0.5694306 0.5514486
+$synth.loglik
+[1] -20.52250
+
> plot(pb)
>
> pb <- probe(
@@ -195,6 +210,9 @@
mean
0.32567433
+$synth.loglik
+[1] -93.0045
+
> plot(pb)
>
> coef(po) <- c(log.r=log(10),log.sigma=log(0.3),log.phi=log(20),N.0=5,e.0=0)
@@ -230,6 +248,9 @@
marg.1 marg.2 marg.3
0.001998002 0.231768232 0.437562438
+$synth.loglik
+[1] -23.38539
+
> plot(pb)
>
> system.time(
@@ -245,15 +266,15 @@
+ )
+ )
user system elapsed
- 10.152 0.147 10.301
+ 13.299 0.053 13.353
> plot(pm)
>
> cbind(truth=coef(ricker),est=coef(pm),guess=coef(po))
truth est guess
-log.r 3.800000 2.635216 2.302585
+log.r 3.800000 2.653141 2.302585
log.sigma -1.203973 -1.203973 -1.203973
-log.phi 2.302585 3.604724 2.995732
-N.0 7.000000 4.967038 5.000000
+log.phi 2.302585 3.676769 2.995732
+N.0 7.000000 3.715260 5.000000
e.0 0.000000 0.000000 0.000000
>
> pb <- probe(
@@ -286,6 +307,9 @@
nlar.1^1 nlar.1^2 nlar.1^3
0.003996004 0.157842158 0.131868132
+$synth.loglik
+[1] -8.830676
+
> plot(pb)
>
> system.time(
@@ -301,17 +325,17 @@
+ )
+ )
user system elapsed
- 6.215 0.093 6.310
+ 8.512 0.041 8.554
> plot(pm)
> plot(as(pm,"pomp"),variables="y")
> plot(simulate(pm),variables="y")
>
> cbind(truth=coef(ricker),est=coef(pm),guess=coef(po))
truth est guess
-log.r 3.800000 1.650922 2.302585
+log.r 3.800000 1.568617 2.302585
log.sigma -1.203973 -1.203973 -1.203973
-log.phi 2.302585 3.424852 2.995732
-N.0 7.000000 5.139255 5.000000
+log.phi 2.302585 3.730399 2.995732
+N.0 7.000000 5.375539 5.000000
e.0 0.000000 0.000000 0.000000
>
> pb <- probe(
@@ -328,7 +352,7 @@
+ )
> pb at datvals
marg.1 marg.2 marg.3
-25.118195 -4.413983 -4.482302
+24.306961 -4.195941 -3.990929
> summary(pb)
$coef
log.r log.sigma log.phi N.0 e.0
@@ -339,12 +363,15 @@
$quantiles
marg.1 marg.2 marg.3
- 1.000 0.000 0.065
+ 1.000 0.000 0.067
$pvals
marg.1 marg.2 marg.3
-0.001998002 0.001998002 0.131868132
+0.001998002 0.001998002 0.135864136
+$synth.loglik
+[1] -66.05758
+
> plot(pb)
>
> pb <- probe(
@@ -377,6 +404,9 @@
acf.0.y acf.1.y acf.2.y acf.3.y acf.4.y acf.5.y
0.8171828 0.1358641 0.6313686 0.1878122 0.2897103 0.6853147
+$synth.loglik
+[1] -14.65429
+
> plot(pb)
>
> pb <- probe(
@@ -409,6 +439,9 @@
acf.1.y acf.2.y acf.3.y acf.4.y acf.5.y
0.3056943 0.6353646 0.2097902 0.2877123 0.6953047
+$synth.loglik
+[1] 5.22575
+
> plot(pb)
>
> pb <- probe(
@@ -450,6 +483,9 @@
v acf.0.y acf.1.y acf.2.y nlar.1^1 nlar.2^1
0.8171828 0.8171828 0.1358641 0.6313686 0.2117882 0.2677323
+$synth.loglik
+[1] NaN
+
> plot(pb)
>
> try(
@@ -473,14 +509,14 @@
+ probes=list(
+ mn=probe.mean("y",transform=sqrt,trim=0.1),
+ md=function(y)median(as.numeric(y)),
-+ wacko=function(y) if (y[1]==69) 1 else c(1,2)
++ wacko=function(y) if (y[1]==68) 1 else c(1,2)
+ ),
+ nsim=100,
+ seed=838775L
+ )
+ )
Error in .local(object, probes, ...) :
- variable-sized results returned by probe 2
+ probes return different number of values on different datasets
>
>
> try(
More information about the pomp-commits
mailing list