[Rcpp-commits] r3512 - in pkg/RcppSMC: . R inst/include man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Mar 17 02:15:56 CET 2012
Author: edd
Date: 2012-03-17 02:15:56 +0100 (Sat, 17 Mar 2012)
New Revision: 3512
Modified:
pkg/RcppSMC/ChangeLog
pkg/RcppSMC/R/pfEx.R
pkg/RcppSMC/TODO
pkg/RcppSMC/inst/include/pffuncs.h
pkg/RcppSMC/man/pfEx.Rd
pkg/RcppSMC/src/blockpfgaussianopt.cpp
pkg/RcppSMC/src/pf.cpp
Log:
reworked pfEx a little to take data from R
Modified: pkg/RcppSMC/ChangeLog
===================================================================
--- pkg/RcppSMC/ChangeLog 2012-03-16 11:09:11 UTC (rev 3511)
+++ pkg/RcppSMC/ChangeLog 2012-03-17 01:15:56 UTC (rev 3512)
@@ -1,4 +1,12 @@
+2012-03-16 Dirk Eddelbuettel <edd at dexter>
+
+ * src/pf.cpp: Changed to get example data from R and pass it to a
+ vector of cv_obs
+ * R/pfEx.R: Added a helper func. to read data, restructured pfEx()
+ * man/pfEx.Rd: Updated accordingly
+
2012-03-16 Adam Johansen <a.m.johansen at warwick.ac.uk>
+
* NAMESPACE: Added pfNonlinBS.
* src/pfnonlinbS.cpp: Bootstrap particle nonlinear particle filter
example.
@@ -18,9 +26,9 @@
* src/pffuncs.cpp merged into src/pf.cpp.
2012-03-13 Dirk Eddelbuettel <edd at debian.org>
-
+
* src/blockpfgaussianoptfuncs.cc: minor new/delete fix
-
+
2012-03-13 Adam Johansen <a.m.johansen at warwick.ac.uk>
* NAMESPACE: Added blockpfGaussianOpt.
Modified: pkg/RcppSMC/R/pfEx.R
===================================================================
--- pkg/RcppSMC/R/pfEx.R 2012-03-16 11:09:11 UTC (rev 3511)
+++ pkg/RcppSMC/R/pfEx.R 2012-03-17 01:15:56 UTC (rev 3512)
@@ -1,21 +1,31 @@
-pfEx<- function(filename="", particles=1000, plot=FALSE) {
+pfEx<- function(data, particles=1000, plot=FALSE) {
- if (filename=="") {
- filename <- system.file("sampleData", "pf-data.csv", package="RcppSMC")
- }
- res <- .Call("pf", filename, particles, package="RcppSMC")
+ # if no data supplied, use default
+ if (missing(data)) data <- getEx1Data()
+ # more eloquent tests can be added
+ stopifnot(nrow(data) > 0,
+ ncol(data) == 2,
+ colnames(data) == c("x", "y"))
+
+ res <- .Call("pf", as.matrix(data), particles, package="RcppSMC")
+
if (plot) {
## plot 5.1 from vignette / paper
- data <- read.table(file=filename, skip=1, header=FALSE,
- col.names=c("x","y"), sep="")
- with(data, plot(x,y,col="red"))
+ with(data, plot(x, y, col="red"))
with(res, lines(Xm, Ym, lty="dashed"))
}
invisible(res)
}
+# simple convenience function, should probably make the data a
+# data component of the package...
+getEx1Data <- function() {
+ file <- system.file("sampleData", "pf-data.csv", package="RcppSMC")
+ dat <- read.table(file, skip=1, header=FALSE, col.names=c("x","y"))
+ invisible(dat)
+}
Modified: pkg/RcppSMC/TODO
===================================================================
--- pkg/RcppSMC/TODO 2012-03-16 11:09:11 UTC (rev 3511)
+++ pkg/RcppSMC/TODO 2012-03-17 01:15:56 UTC (rev 3512)
@@ -3,7 +3,10 @@
* Maybe remove example 2 from the JSS paper
- * For blockpfGaussianOpt, we could be fancy and return an S3 object
+ * Maybe rename example 1 from the JSS paper -- what would be a good
+ name instead of 'pfEx'? Idem for the C++ code.
+
+ * For blockpfGaussianOpt et al, we could be fancy and return an S3 object
with a plot method. Maybe later.
* Once we are happy with the GSL -> R rng switch, we can remove all
Modified: pkg/RcppSMC/inst/include/pffuncs.h
===================================================================
--- pkg/RcppSMC/inst/include/pffuncs.h 2012-03-16 11:09:11 UTC (rev 3511)
+++ pkg/RcppSMC/inst/include/pffuncs.h 2012-03-17 01:15:56 UTC (rev 3512)
@@ -25,4 +25,6 @@
extern double nu_y;
extern double Delta;
-extern cv_obs * y;
+//extern cv_obs * y;
+//extern Rcpp::NumericMatrix y;
+extern std::vector<cv_obs> y;
Modified: pkg/RcppSMC/man/pfEx.Rd
===================================================================
--- pkg/RcppSMC/man/pfEx.Rd 2012-03-16 11:09:11 UTC (rev 3511)
+++ pkg/RcppSMC/man/pfEx.Rd 2012-03-17 01:15:56 UTC (rev 3512)
@@ -8,13 +8,12 @@
tracking' problem of 100 observations is solved with 1000 particles.
}
\usage{
- pfEx(filename="", particles=1000, plot=FALSE)
+ pfEx(data, particles=1000, plot=FALSE)
}
\arguments{
- \item{filename}{A character variable describing path and filename of a
- data file in the format described in the \pkg{SMCTC}
- documentation. Note we may replace this with a direct transfer of a
- data matrix.}
+ \item{data}{A two-column matrix or dataframe containing x and y
+ values. The default data set from Johansen (2009) is used as the
+ default if no data is supplied.}
\item{particles}{An integer specifying the number of particles.}
\item{plot}{A boolean variable describing whether plot should
illustrate the estimated path along with the data.}
Modified: pkg/RcppSMC/src/blockpfgaussianopt.cpp
===================================================================
--- pkg/RcppSMC/src/blockpfgaussianopt.cpp 2012-03-16 11:09:11 UTC (rev 3511)
+++ pkg/RcppSMC/src/blockpfgaussianopt.cpp 2012-03-17 01:15:56 UTC (rev 3512)
@@ -11,7 +11,7 @@
///The observations
namespace BSPFG {
Rcpp::NumericVector y;
-};
+}
long lLag = 1;
Modified: pkg/RcppSMC/src/pf.cpp
===================================================================
--- pkg/RcppSMC/src/pf.cpp 2012-03-16 11:09:11 UTC (rev 3511)
+++ pkg/RcppSMC/src/pf.cpp 2012-03-17 01:15:56 UTC (rev 3512)
@@ -39,8 +39,10 @@
using namespace std;
///The observations
-cv_obs * y;
-long load_data(char const * szName, cv_obs** y);
+//cv_obs * y;
+//Rcpp::NumericMatrix y;
+std::vector<cv_obs> y;
+//long load_data(char const * szName, cv_obs** y);
double integrand_mean_x(const cv_state&, void*);
double integrand_mean_y(const cv_state&, void*);
@@ -48,17 +50,25 @@
double integrand_var_y(const cv_state&, void*);
// pf() function callable from R via Rcpp:: essentially the same as main() from pf.cc
-extern "C" SEXP pf(SEXP fileS, SEXP partS) {
+// minor interface change to pass data down as matrix, rather than a filename
+extern "C" SEXP pf(SEXP dataS, SEXP partS) {
long lIterates;
try {
- std::string filename = Rcpp::as<std::string>(fileS);
+ //std::string filename = Rcpp::as<std::string>(fileS);
unsigned long lNumber = Rcpp::as<unsigned long>(partS);
- //Load observations
- lIterates = load_data(filename.c_str(), &y);
+ // Load observations -- or rather copy them in from R
+ //lIterates = load_data(filename.c_str(), &y);
+ Rcpp::NumericMatrix dat = Rcpp::NumericMatrix(dataS); // so we expect a matrix
+ lIterates = dat.nrow();
+ y.reserve(lIterates);
+ for (long i = 0; i < lIterates; ++i) {
+ y[i].x_pos = dat(i,0);
+ y[i].y_pos = dat(i,1);
+ }
//Initialise and run the sampler
smc::sampler<cv_state> Sampler(lNumber, SMC_HISTORY_NONE);
@@ -91,6 +101,7 @@
return R_NilValue; // to provide a return
}
+#if 0
long load_data(char const * szName, cv_obs** yp)
{
FILE * fObs = fopen(szName,"rt");
@@ -114,6 +125,7 @@
return lIterates;
}
+#endif
double integrand_mean_x(const cv_state& s, void *)
{
@@ -160,7 +172,13 @@
/// \param X The state to consider
double logLikelihood(long lTime, const cv_state & X)
{
- return - 0.5 * (nu_y + 1.0) * (log(1 + pow((X.x_pos - y[lTime].x_pos)/scale_y,2) / nu_y) + log(1 + pow((X.y_pos - y[lTime].y_pos)/scale_y,2) / nu_y));
+ //return - 0.5 * (nu_y + 1.0) * (log(1 + pow((X.x_pos - y[lTime].x_pos)/scale_y,2) / nu_y) + log(1 + pow((X.y_pos - y[lTime].y_pos)/scale_y,2) / nu_y));
+
+ //double y_xpos = y(lTime,0);
+ //double y_ypos = y(lTime,1);
+ //return - 0.5 * (nu_y + 1.0) * (log(1 + pow((X.x_pos - y_xpos)/scale_y,2) / nu_y) + log(1 + pow((X.y_pos - y_ypos)/scale_y,2) / nu_y));
+
+ return - 0.5 * (nu_y + 1.0) * (log(1 + pow((X.x_pos - y[lTime].x_pos)/scale_y,2) / nu_y) + log(1 + pow((X.y_pos - y[lTime].y_pos)/scale_y,2) / nu_y));
}
///A function to initialise particles
More information about the Rcpp-commits
mailing list