[Pomp-commits] r450 - in pkg: . R inst inst/doc src tests
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue May 3 17:41:13 CEST 2011
Author: kingaa
Date: 2011-05-03 17:41:13 +0200 (Tue, 03 May 2011)
New Revision: 450
Removed:
pkg/src/resample.c
Modified:
pkg/DESCRIPTION
pkg/R/pfilter.R
pkg/inst/ChangeLog
pkg/inst/doc/advanced_topics_in_pomp.pdf
pkg/inst/doc/intro_to_pomp.Rnw
pkg/inst/doc/intro_to_pomp.pdf
pkg/src/pfilter.c
pkg/tests/ou2-probe.Rout.save
pkg/tests/ricker-probe.Rout.save
pkg/tests/ricker.Rout.save
Log:
- move resampling and random-walking inside 'pfilter.internal' into C (pfilter.c)
- move resampling codes from 'resample.c' to 'pfilter.c'
- update test .save files
- minor corrections in 'intro_to_pomp' vignette
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2011-05-03 14:52:19 UTC (rev 449)
+++ pkg/DESCRIPTION 2011-05-03 15:41:13 UTC (rev 450)
@@ -1,8 +1,8 @@
Package: pomp
Type: Package
Title: Statistical inference for partially observed Markov processes
-Version: 0.36-5
-Date: 2011-04-24
+Version: 0.36-6
+Date: 2011-05-04
Author: Aaron A. King, Edward L. Ionides, Carles Breto, Steve Ellner, Bruce Kendall, Helen Wearing,
Matthew J. Ferrari, Michael Lavine, Daniel C. Reuman
Maintainer: Aaron A. King <kingaa at umich.edu>
Modified: pkg/R/pfilter.R
===================================================================
--- pkg/R/pfilter.R 2011-05-03 14:52:19 UTC (rev 449)
+++ pkg/R/pfilter.R 2011-05-03 15:41:13 UTC (rev 450)
@@ -87,6 +87,7 @@
sigma <- .rw.sd
} else {
rw.names <- character(0)
+ sigma <- NULL
}
loglik <- rep(NA,ntimes)
@@ -189,23 +190,27 @@
}
## prediction mean, prediction variance, filtering mean, effective sample size, log-likelihood
+ ## also do resampling if filtering has not failed
xx <- try(
.Call(
pfilter_computations,
X,params,
- random.walk,rw.names,
+ random.walk,sigma,
pred.mean,pred.var,
filter.mean,weights,tol
),
silent=FALSE
)
if (inherits(xx,'try-error')) {
- stop(sQuote("pfilter")," error: error in prediction mean/variance computation",call.=FALSE)
+ stop(sQuote("pfilter")," error",call.=FALSE)
}
all.fail <- xx$fail
loglik[nt] <- xx$loglik
eff.sample.size[nt] <- xx$ess
+ x <- xx$states
+ params <- xx$params
+
if (pred.mean)
pred.m[,nt] <- xx$pm
if (pred.var)
@@ -219,19 +224,8 @@
message("filtering failure at time t = ",times[nt+1])
if (nfail>max.fail)
stop(sQuote("pfilter")," error: too many filtering failures",call.=FALSE)
- x[,] <- X[,,1,drop=FALSE]
- } else { ## matrix with samples (columns) from filtering distribution theta.t | Y.t
- sample <- .Call(systematic_resampling,weights)
- x[,] <- X[,sample,1,drop=FALSE]
- params <- params[,sample,drop=FALSE]
}
-
- ## random walk for parameters
- if (random.walk) {
- pred.v[rw.names,nt] <- pred.v[rw.names,nt]+sigma^2
- params[rw.names,] <- rnorm(n=Np*length(sigma),mean=params[rw.names,],sd=sigma)
- }
-
+
if (save.states) {
xparticles[,,nt] <- x
}
Modified: pkg/inst/ChangeLog
===================================================================
--- pkg/inst/ChangeLog 2011-05-03 14:52:19 UTC (rev 449)
+++ pkg/inst/ChangeLog 2011-05-03 15:41:13 UTC (rev 450)
@@ -1,3 +1,12 @@
+2011-05-03 kingaa
+
+ * [r449] R/bsplines.R, man/bsplines.Rd: - add 'names' argument
+
+2011-04-27 kingaa
+
+ * [r448] inst/ChangeLog, inst/NEWS: - version 0.36-5
+ - update NEWS file and Changelog
+
2011-04-24 kingaa
* [r447] inst/ChangeLog, inst/doc/advanced_topics_in_pomp.pdf,
Modified: pkg/inst/doc/advanced_topics_in_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/inst/doc/intro_to_pomp.Rnw
===================================================================
--- pkg/inst/doc/intro_to_pomp.Rnw 2011-05-03 14:52:19 UTC (rev 449)
+++ pkg/inst/doc/intro_to_pomp.Rnw 2011-05-03 15:41:13 UTC (rev 450)
@@ -140,14 +140,14 @@
Later we'll look at a continuous-time model for which no such special tricks are available.
Consider the discrete-time Gompertz model of population growth.
-Under this model, the density, $X_{t+1}$, of a population of plants or animals at time $t+{\Delta}t$ depends on the density, $X_{t}$, at time $t$ according to
+Under this model, the density, $X_{t+{\Delta}t}$, of a population of plants or animals at time $t+{\Delta}t$ depends on the density, $X_{t}$, at time $t$ according to
\begin{equation}\label{eq:gompertz1}
- X_{t+1}=K^{1-e^{-r\,{\Delta}t}}\,X_{t}^{e^{-r\,{\Delta}t}}\,\varepsilon_{t},
+ X_{t+{\Delta}t}=K^{1-e^{-r\,{\Delta}t}}\,X_{t}^{e^{-r\,{\Delta}t}}\,\varepsilon_{t},
\end{equation}
where $K$ is the so-called ``carrying capacity'' of the population, $r$ is a positive parameter, and the $\varepsilon_{t}$ are independent and identically-distributed lognormal random variables.
In different notation, this model is
\begin{equation}\label{eq:gompertz2}
- \log{X_{t+1}}\;\sim\;\mathrm{normal}(\log{K}+\log{\left(\frac{X_t}{K}\right)}\,e^{-r\,{\Delta}t},\sigma),
+ \log{X_{t+{\Delta}t}}\;\sim\;\mathrm{normal}(\log{K}+\log{\left(\frac{X_t}{K}\right)}\,e^{-r\,{\Delta}t},\sigma),
\end{equation}
where $\sigma^2={\mathrm{Var}[\log{\epsilon_{t}}]}$.
We'll assume that we can measure the population density only with error.
@@ -166,12 +166,12 @@
In order to fully specify this partially-observed Markov process, we must implement both the process model (i.e., the unobserved process) and the measurement model (the observation process).
As we saw before, we would like to be able to:
\begin{enumerate}
-\item \label{it:rproc} simulate from the process model, i.e., make a random draw from $X_{t+1}\,\vert\,X_{t}=x$ for arbitrary $x$ and $t$ (\code{rprocess}),
-\item \label{it:dproc} compute the probability density function (pdf) of state transitions, i.e., compute $f(X_{t+1}=x'\,\vert\,X_{t}=x)$ for arbitrary $x$, $x'$, and $t$ (\code{dprocess}),
+\item \label{it:rproc} simulate from the process model, i.e., make a random draw from $X_{t+{\Delta}t}\,\vert\,X_{t}=x$ for arbitrary $x$ and $t$ (\code{rprocess}),
+\item \label{it:dproc} compute the probability density function (pdf) of state transitions, i.e., compute $f(X_{t+{\Delta}t}=x'\,\vert\,X_{t}=x)$ for arbitrary $x$, $x'$, $t$, and ${\Delta}t$ (\code{dprocess}),
\item \label{it:rmeas} simulate from the measurement model, i.e., make a random draw from $Y_{t}\,\vert\,X_{t}=x$ for arbitrary $x$ and $t$ (\code{rmeasure}),
\item \label{it:dmeas} compute the measurement model pdf, i.e., $f(Y_{t}=y\,\vert\,X_{t}=x)$ for arbitrary $x$, $y$, and $t$ (\code{dmeasure}), and
\item \label{it:skel} compute the \emph{deterministic skeleton}.
- In discrete-time, this is the map $x\,\mapsto\,\mathbb{E}[X_{t+1}\,\vert\,X_{t}=x]$ for arbitrary $x$.
+ In discrete-time, this is the map $x\,\mapsto\,\mathbb{E}[X_{t+{\Delta}t}\,\vert\,X_{t}=x]$ for arbitrary $x$.
\end{enumerate}
For this simple model, all this is easy enough.
More generally, it will be difficult to do some of these things.
@@ -257,7 +257,7 @@
Note that \code{gompertz.meas.dens} is closely related to \code{gompertz.meas.sim}.
\clearpage
-\section{Simulating the model}
+\section{Simulating the model: \code{simulate}}
With the two functions above, we already have all we need to simulate the full model.
The first step is to construct an \R\ object of class \pomp\ which will serve as a container to hold the model and data.
@@ -342,7 +342,7 @@
\label{fig:gompertz-first-simulation-plot}
\end{figure}
-\section{Computing likelihood using particle filtering}
+\section{Computing likelihood using particle filtering: \code{pfilter}}
Some parameter estimation algorithms in the \pomp\ package only require \code{rprocess} and \code{rmeasure}.
These include the nonlinear forecasting algorithm \code{nlf} and the probe-matching algorithm \code{probe.match}.
Modified: pkg/inst/doc/intro_to_pomp.pdf
===================================================================
(Binary files differ)
Modified: pkg/src/pfilter.c
===================================================================
--- pkg/src/pfilter.c 2011-05-03 14:52:19 UTC (rev 449)
+++ pkg/src/pfilter.c 2011-05-03 15:41:13 UTC (rev 450)
@@ -3,6 +3,45 @@
#include "pomp_internal.h"
#include <Rdefines.h>
+static void nosort_resamp (int n, double *w, int *p, int offset)
+{
+ int i, j;
+ double du, u;
+
+ for (j = 1; j < n; j++) w[j] += w[j-1];
+
+ if (w[n-1] <= 0.0)
+ error("non-positive sum of weights");
+
+ du = w[n-1] / ((double) n);
+ u = runif(-du,0);
+
+ for (i = 0, j = 0; j < n; j++) {
+ u += du;
+ while (u > w[i]) i++;
+ p[j] = i;
+ }
+ if (offset) // add offset if needed
+ for (j = 0; j < n; j++) p[j] += offset;
+
+}
+
+SEXP systematic_resampling (SEXP weights)
+{
+ int nprotect = 0;
+ double u, du, *w;
+ int i, j, *p, n;
+ SEXP perm;
+
+ GetRNGstate();
+ n = LENGTH(weights);
+ PROTECT(perm = NEW_INTEGER(n)); nprotect++;
+ nosort_resamp(n,REAL(weights),INTEGER(perm),1);
+ UNPROTECT(nprotect);
+ PutRNGstate();
+ return(perm);
+}
+
// examines weights for filtering failure
// computes log likelihood and effective sample size
// computes (if desired) prediction mean, prediction variance, filtering mean.
@@ -11,15 +50,17 @@
// if length(weights) == 1, an unweighted average is computed.
// returns all of the above in a named list
SEXP pfilter_computations (SEXP x, SEXP params,
- SEXP rw, SEXP rw_names,
+ SEXP rw, SEXP rw_sd,
SEXP predmean, SEXP predvar,
SEXP filtmean, SEXP weights, SEXP tol)
{
int nprotect = 0;
- SEXP pm = R_NilValue, pv = R_NilValue, fm = R_NilValue, ess, fail, loglik;
+ SEXP pm = R_NilValue, pv = R_NilValue, fm = R_NilValue;
+ SEXP rw_names, ess, fail, loglik;
+ SEXP newstates, newparams;
SEXP retval, retvalnames;
double *xpm = 0, *xpv = 0, *xfm = 0, *xw = 0, *xx = 0, *xp = 0;
- SEXP dimX, Pnames, pindex;
+ SEXP dimX, dimP, newdim, Xnames, Pnames, pindex;
int *dim, *pidx, lv;
int nvars, npars = 0, nrw = 0, nreps, offset, nlost;
int do_rw, do_pm, do_pv, do_fm, all_fail = 0;
@@ -30,18 +71,29 @@
dim = INTEGER(dimX);
nvars = dim[0]; nreps = dim[1];
xx = REAL(x);
- do_rw = *(LOGICAL(AS_LOGICAL(rw)));
- do_pm = *(LOGICAL(AS_LOGICAL(predmean)));
- do_pv = *(LOGICAL(AS_LOGICAL(predvar)));
- do_fm = *(LOGICAL(AS_LOGICAL(filtmean)));
+ PROTECT(Xnames = GET_ROWNAMES(GET_DIMNAMES(x))); nprotect++;
- PROTECT(ess = NEW_NUMERIC(1)); nprotect++;
- PROTECT(loglik = NEW_NUMERIC(1)); nprotect++;
- PROTECT(fail = NEW_LOGICAL(1)); nprotect++;
+ PROTECT(dimP = GET_DIM(params)); nprotect++;
+ dim = INTEGER(dimP);
+ npars = dim[0];
+ if (nreps != dim[1])
+ error("'states' and 'params' do not agree in second dimension");
+ PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
- xw = REAL(weights);
- toler = *(REAL(tol));
+ PROTECT(rw_names = GET_NAMES(rw_sd)); nprotect++; // names of parameters undergoing random walk
+ do_rw = *(LOGICAL(AS_LOGICAL(rw))); // do random walk in parameters?
+ do_pm = *(LOGICAL(AS_LOGICAL(predmean))); // calculate prediction means?
+ do_pv = *(LOGICAL(AS_LOGICAL(predvar))); // calculate prediction variances?
+ do_fm = *(LOGICAL(AS_LOGICAL(filtmean))); // calculate filtering means?
+
+ PROTECT(ess = NEW_NUMERIC(1)); nprotect++; // effective sample size
+ PROTECT(loglik = NEW_NUMERIC(1)); nprotect++; // log likelihood
+ PROTECT(fail = NEW_LOGICAL(1)); nprotect++; // particle failure?
+
+ xw = REAL(weights);
+ toler = *(REAL(tol)); // failure tolerance
+
// check the weights and compute sum and sum of squares
for (k = 0, w = 0, ws = 0, nlost = 0; k < nreps; k++) {
if (xw[k] > toler) {
@@ -55,19 +107,17 @@
if (nlost >= nreps) all_fail = 1; // all particles are lost
if (all_fail) {
*(REAL(loglik)) = log(toler); // minimum log-likelihood
- *(REAL(ess)) = 0; // zero effective sample size
+ *(REAL(ess)) = 0; // zero effective sample size
} else {
*(REAL(loglik)) = log(w/((double) nreps)); // mean of weights is likelihood
- *(REAL(ess)) = w*w/ws; // effective sample size
+ *(REAL(ess)) = w*w/ws; // effective sample size
}
*(LOGICAL(fail)) = all_fail;
if (do_rw) {
- PROTECT(Pnames = GET_ROWNAMES(GET_DIMNAMES(params))); nprotect++;
- PROTECT(pindex = matchnames(Pnames,rw_names)); nprotect++;
+ PROTECT(pindex = matchnames(Pnames,rw_names)); nprotect++; // indices of parameters undergoing random walk
xp = REAL(params);
pidx = INTEGER(pindex);
- npars = LENGTH(Pnames);
nrw = LENGTH(rw_names);
lv = nvars+nrw;
} else {
@@ -114,10 +164,10 @@
// compute filter mean
if (do_fm) {
- if (all_fail) {
+ if (all_fail) { // unweighted average
for (k = 0, ws = 0; k < nreps; k++) ws += xx[j+k*nvars];
xfm[j] = ws/((double) nreps);
- } else {
+ } else { // weighted average
for (k = 0, ws = 0; k < nreps; k++) ws += xx[j+k*nvars]*xw[k];
xfm[j] = ws/w;
}
@@ -125,10 +175,10 @@
}
- // compute means and variances for parameters
+ // compute means and variances for parameters (if needed)
if (do_rw) {
for (j = 0; j < nrw; j++) {
- offset = pidx[j];
+ offset = pidx[j]; // position of the parameter
if (do_pm || do_pv) {
for (k = 0, sum = 0; k < nreps; k++) sum += xp[offset+k*npars];
@@ -145,12 +195,13 @@
}
}
+
if (do_fm) {
for (j = 0; j < npars; j++) {
- if (all_fail) {
+ if (all_fail) { // unweighted average
for (k = 0, ws = 0; k < nreps; k++) ws += xp[j+k*npars];
xfm[nvars+j] = ws/((double) nreps);
- } else {
+ } else { // weighted average
for (k = 0, ws = 0; k < nreps; k++) ws += xp[j+k*npars]*xw[k];
xfm[nvars+j] = ws/w;
}
@@ -158,28 +209,89 @@
}
}
- PROTECT(retval = NEW_LIST(6)); nprotect++;
- PROTECT(retvalnames = NEW_CHARACTER(6)); nprotect++;
+ GetRNGstate();
+
+ if (!all_fail) { // resample the particles unless we have filtering failure
+ int xdim[2];
+ int sample[nreps];
+ double *ss, *st, *ps, *pt;
+
+ // create storage for new states
+ xdim[0] = nvars; xdim[1] = nreps;
+ PROTECT(newstates = makearray(2,xdim)); nprotect++;
+ setrownames(newstates,Xnames,2);
+ ss = REAL(x);
+ st = REAL(newstates);
+
+ // create storage for new parameters
+ xdim[0] = npars; xdim[1] = nreps;
+ PROTECT(newparams = makearray(2,xdim)); nprotect++;
+ setrownames(newparams,Pnames,2);
+ ps = REAL(params);
+ pt = REAL(newparams);
+
+ // resample
+ nosort_resamp(nreps,REAL(weights),sample,0);
+ for (k = 0; k < nreps; k++) { // copy the particles
+ for (j = 0, xx = ss+nvars*sample[k]; j < nvars; j++, st++, xx++) *st = *xx;
+ for (j = 0, xp = ps+npars*sample[k]; j < npars; j++, pt++, xp++) *pt = *xp;
+ }
+ } else { // don't resample: just drop 3rd dimension in x prior to return
+ PROTECT(newdim = NEW_INTEGER(2)); nprotect++;
+ dim = INTEGER(newdim);
+ dim[0] = nvars; dim[1] = nreps;
+ SET_DIM(x,newdim);
+ setrownames(x,Xnames,2);
+ }
+
+ if (do_rw) { // if random walk, adjust prediction variance and move particles
+ xx = REAL(rw_sd);
+ xp = (all_fail) ? REAL(params) : REAL(newparams);
+ for (j = 0; j < nrw; j++) {
+ offset = pidx[j];
+ vsq = xx[j];
+ if (do_pv) {
+ xpv[nvars+j] += vsq*vsq;
+ }
+ for (k = 0; k < nreps; k++)
+ xp[offset+k*npars] += rnorm(0,vsq);
+ }
+ }
+
+ PutRNGstate();
+
+ PROTECT(retval = NEW_LIST(8)); nprotect++;
+ PROTECT(retvalnames = NEW_CHARACTER(8)); nprotect++;
SET_STRING_ELT(retvalnames,0,mkChar("fail"));
SET_STRING_ELT(retvalnames,1,mkChar("loglik"));
SET_STRING_ELT(retvalnames,2,mkChar("ess"));
- SET_STRING_ELT(retvalnames,3,mkChar("pm"));
- SET_STRING_ELT(retvalnames,4,mkChar("pv"));
- SET_STRING_ELT(retvalnames,5,mkChar("fm"));
+ SET_STRING_ELT(retvalnames,3,mkChar("states"));
+ SET_STRING_ELT(retvalnames,4,mkChar("params"));
+ SET_STRING_ELT(retvalnames,5,mkChar("pm"));
+ SET_STRING_ELT(retvalnames,6,mkChar("pv"));
+ SET_STRING_ELT(retvalnames,7,mkChar("fm"));
SET_NAMES(retval,retvalnames);
SET_ELEMENT(retval,0,fail);
SET_ELEMENT(retval,1,loglik);
SET_ELEMENT(retval,2,ess);
+
+ if (all_fail) {
+ SET_ELEMENT(retval,3,x);
+ SET_ELEMENT(retval,4,params);
+ } else {
+ SET_ELEMENT(retval,3,newstates);
+ SET_ELEMENT(retval,4,newparams);
+ }
if (do_pm) {
- SET_ELEMENT(retval,3,pm);
+ SET_ELEMENT(retval,5,pm);
}
if (do_pv) {
- SET_ELEMENT(retval,4,pv);
+ SET_ELEMENT(retval,6,pv);
}
if (do_fm) {
- SET_ELEMENT(retval,5,fm);
+ SET_ELEMENT(retval,7,fm);
}
UNPROTECT(nprotect);
Deleted: pkg/src/resample.c
===================================================================
--- pkg/src/resample.c 2011-05-03 14:52:19 UTC (rev 449)
+++ pkg/src/resample.c 2011-05-03 15:41:13 UTC (rev 450)
@@ -1,93 +0,0 @@
-// -*- C++ -*-
-
-#include "pomp_internal.h"
-#include <Rdefines.h>
-
-static void nosort_resamp (int n, double *w, int *p, int offset)
-{
- int i, j;
- double du, u;
-
- for (j = 1; j < n; j++) w[j] += w[j-1];
-
- if (w[n-1] <= 0.0)
- error("non-positive sum of weights");
-
- GetRNGstate();
- du = w[n-1] / ((double) n);
- u = runif(-du,0);
- PutRNGstate();
-
- for (i = 0, j = 0; j < n; j++) {
- u += du;
- while (u > w[i]) i++;
- p[j] = i;
- }
- if (offset) // add offset if needed
- for (j = 0; j < n; j++) p[j] += offset;
-
-}
-
-SEXP systematic_resampling (SEXP weights)
-{
- int nprotect = 0;
- double u, du, *w;
- int i, j, *p, n;
- SEXP perm;
-
- n = LENGTH(weights);
- PROTECT(perm = NEW_INTEGER(n)); nprotect++;
- nosort_resamp(n,REAL(weights),INTEGER(perm),1);
- UNPROTECT(nprotect);
- return(perm);
-}
-
-SEXP pfilter_resample (SEXP weights, SEXP states, SEXP params)
-{
- int nprotect = 0;
- int xdim[2], nvars, npars, nreps;
- SEXP sdim, pdim, newstates, newparams, sample;
- SEXP retval, retvalnames;
- int i, j, n;
- double *ss, *st, *ps, *pt, *xx;
-
- PROTECT(sdim = GET_DIM(states)); nprotect++;
- nvars = INTEGER(sdim)[0]; nreps = INTEGER(sdim)[1];
-
- xdim[0] = nvars; xdim[1] = nreps;
- PROTECT(newstates = makearray(2,xdim)); nprotect++;
- setrownames(newstates,GET_ROWNAMES(GET_DIMNAMES(states)),2);
- ss = REAL(states);
- st = REAL(newstates);
-
- PROTECT(pdim = GET_DIM(params)); nprotect++;
- npars = INTEGER(pdim)[0];
- if (nreps != INTEGER(pdim)[1])
- error("'states' and 'params' do not agree in second dimension");
-
- xdim[0] = npars; xdim[1] = nreps;
- PROTECT(newparams = makearray(2,xdim)); nprotect++;
- setrownames(newparams,GET_ROWNAMES(GET_DIMNAMES(params)),2);
- ps = REAL(params);
- pt = REAL(newparams);
-
- {
- int sample[nreps];
- nosort_resamp(nreps,REAL(weights),sample,0);
- for (i = 0; i < nreps; i++) {
- for (j = 0, xx = ss+nvars*sample[i]; j < nvars; j++, st++, xx++) *st = *xx;
- for (j = 0, xx = ps+npars*sample[i]; j < npars; j++, pt++, xx++) *pt = *xx;
- }
- }
-
- PROTECT(retval = NEW_LIST(2)); nprotect++;
- SET_ELEMENT(retval,0,newstates);
- SET_ELEMENT(retval,1,newparams);
- PROTECT(retvalnames = NEW_CHARACTER(2)); nprotect++;
- SET_STRING_ELT(retvalnames,0,mkChar("states"));
- SET_STRING_ELT(retvalnames,1,mkChar("params"));
- SET_NAMES(retval,retvalnames);
-
- UNPROTECT(nprotect);
- return retval;
-}
Modified: pkg/tests/ou2-probe.Rout.save
===================================================================
--- pkg/tests/ou2-probe.Rout.save 2011-05-03 14:52:19 UTC (rev 449)
+++ pkg/tests/ou2-probe.Rout.save 2011-05-03 15:41:13 UTC (rev 450)
@@ -1,7 +1,8 @@
-R version 2.11.1 (2010-05-31)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.12.2 (2011-02-25)
+Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
+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.
@@ -16,6 +17,9 @@
Type 'q()' to quit R.
> library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
> set.seed(1066L)
>
> data(ou2)
@@ -189,7 +193,7 @@
> y2 <- head(x$y2,-1)
> z2 <- tail(x$y2,-1)
> max(abs(pb at datvals-c(mean(y1*z1)/mean(x$y1^2),mean(y2*z2)/mean(x$y2^2),mean(y1*z1)/mean(y1*y1),mean(y2*z2)/mean(y2*y2))))
-[1] 4.440892e-16
+[1] 1.221245e-15
>
> po <- simulate(ou2)
> pb <- probe(
@@ -219,7 +223,7 @@
0.3016983 0.2957043 0.3216783 0.3016983 0.2957043 0.3216783
$synth.loglik
-[1] 86.8972
+[1] 89.34863
>
> pb <- probe(
Modified: pkg/tests/ricker-probe.Rout.save
===================================================================
--- pkg/tests/ricker-probe.Rout.save 2011-05-03 14:52:19 UTC (rev 449)
+++ pkg/tests/ricker-probe.Rout.save 2011-05-03 15:41:13 UTC (rev 450)
@@ -1,7 +1,8 @@
-R version 2.11.1 (2010-05-31)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.12.2 (2011-02-25)
+Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
+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.
@@ -16,6 +17,9 @@
Type 'q()' to quit R.
> library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
>
> data(ricker)
>
@@ -308,7 +312,7 @@
+ )
+ )
user system elapsed
- 13.482 0.016 13.516
+ 13.888 0.047 13.965
> plot(pm)
>
> cbind(truth=coef(ricker),est=coef(pm),guess=coef(po))
@@ -367,7 +371,7 @@
+ )
+ )
user system elapsed
- 10.165 0.030 10.208
+ 10.304 0.041 10.366
> plot(pm)
> plot(as(pm,"pomp"),variables="y")
> plot(simulate(pm),variables="y")
@@ -526,7 +530,7 @@
0.8171828 0.8171828 0.1358641 0.6313686 0.2117882 0.2677323
$synth.loglik
-[1] 29.83407
+[1] 30.22958
> plot(pb)
>
Modified: pkg/tests/ricker.Rout.save
===================================================================
--- pkg/tests/ricker.Rout.save 2011-05-03 14:52:19 UTC (rev 449)
+++ pkg/tests/ricker.Rout.save 2011-05-03 15:41:13 UTC (rev 450)
@@ -1,7 +1,8 @@
-R version 2.11.1 (2010-05-31)
-Copyright (C) 2010 The R Foundation for Statistical Computing
+R version 2.12.2 (2011-02-25)
+Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
+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.
@@ -16,6 +17,9 @@
Type 'q()' to quit R.
> library(pomp)
+Loading required package: mvtnorm
+Loading required package: subplex
+Loading required package: deSolve
>
> data(ricker)
>
@@ -251,13 +255,13 @@
>
> print(res,digits=3)
guess truth MLE PM
-log.r 3.00 3.80 3.81 3.74
-log.sigma 0.00 -1.20 -1.65 -1.12
+log.r 3.00 3.80 3.79 3.74
+log.sigma 0.00 -1.20 -1.74 -1.12
log.phi 3.00 2.30 2.32 2.41
N.0 7.00 7.00 7.00 7.00
e.0 0.00 0.00 0.00 0.00
-loglik -230.86 -139.55 -137.17 -143.33
-synth.loglik -12.28 17.47 17.59 20.34
+loglik -230.86 -139.55 -137.35 -143.33
+synth.loglik -12.28 17.47 16.47 20.34
>
> plot(ricker)
> plot(simulate(guess))
More information about the pomp-commits
mailing list