[Rcpp-commits] r3525 - in pkg/RcppSMC: . R inst/include man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Mar 18 21:37:55 CET 2012
Author: adamj
Date: 2012-03-18 21:37:55 +0100 (Sun, 18 Mar 2012)
New Revision: 3525
Added:
pkg/RcppSMC/R/pfLineartBS.R
pkg/RcppSMC/inst/include/pflineart.h
pkg/RcppSMC/man/pfLineartBS.Rd
pkg/RcppSMC/src/pflineart.cpp
Removed:
pkg/RcppSMC/R/pfEx.R
pkg/RcppSMC/inst/include/pffuncs.h
pkg/RcppSMC/man/pfEx.Rd
pkg/RcppSMC/src/pf.cpp
Modified:
pkg/RcppSMC/ChangeLog
pkg/RcppSMC/NAMESPACE
Log:
Standardising naming convention for examples.
Modified: pkg/RcppSMC/ChangeLog
===================================================================
--- pkg/RcppSMC/ChangeLog 2012-03-18 19:52:08 UTC (rev 3524)
+++ pkg/RcppSMC/ChangeLog 2012-03-18 20:37:55 UTC (rev 3525)
@@ -1,6 +1,12 @@
2012-03-18 Adam Johansen <a.m.johansen at warwick.ac.uk>
- * src/bspfGaussianOptimal.R tweaked & standardised plotting.
- * src/pfNonlinBS.R fixed plotting bug.
+ * src/pfEx.cpp renamed to pflineartbs.cpp
+ * inst/include/pffuncs.h renamed to pflineartbs.h
+ * R/pfEx.R renamed to pfLineartBS.R
+ * man/pfEx.Rd renamed to pfLineartBS.Rd
+
+ * R/bspfGaussianOptimal.R tweaked & standardised plotting.
+ * R/pfNonlinBS.R fixed plotting bug.
+
* deprecated moved all rareEvents files here.
2012-03-17 Adam Johansen <a.m.johansen at warwick.ac.uk>
Modified: pkg/RcppSMC/NAMESPACE
===================================================================
--- pkg/RcppSMC/NAMESPACE 2012-03-18 19:52:08 UTC (rev 3524)
+++ pkg/RcppSMC/NAMESPACE 2012-03-18 20:37:55 UTC (rev 3525)
@@ -2,7 +2,7 @@
useDynLib(RcppSMC)
export(blockpfGaussianOpt,
- pfEx,
+ pfLineartBS,
pfExOnlinePlot,
pfNonlinBS)
Deleted: pkg/RcppSMC/R/pfEx.R
===================================================================
--- pkg/RcppSMC/R/pfEx.R 2012-03-18 19:52:08 UTC (rev 3524)
+++ pkg/RcppSMC/R/pfEx.R 2012-03-18 20:37:55 UTC (rev 3525)
@@ -1,50 +0,0 @@
-
-pfEx<- function(data, particles=1000, plot=FALSE, onlinePlot) {
-
- # if no data supplied, use default
- if (missing(data)) data <- getEx1Data()
-
- if (missing(onlinePlot)) {
- useOnline <- FALSE
- onlinePlot <- function() { NULL }
- } else {
- useOnline <- TRUE
- # set up x11
- x11(width=3,height=3)
- par(mar=c(3,3,1,1),cex=0.8, pch=19)
- }
-
- # more eloquent tests can be added
- stopifnot(nrow(data) > 0,
- ncol(data) == 2,
- colnames(data) == c("x", "y"),
- class(onlinePlot) == "function")
-
- res <- .Call("pf", as.matrix(data),
- particles,
- useOnline,
- onlinePlot,
- package="RcppSMC")
-
- if (plot) {
- ## plot 5.1 from vignette / paper
- 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)
-}
-
-pfExOnlinePlot <- function(xm, ym) {
- plot(xm, ym, xlim=c(-7,0), ylim=c(2,14))
- Sys.sleep(0.05)
-}
-
Copied: pkg/RcppSMC/R/pfLineartBS.R (from rev 3519, pkg/RcppSMC/R/pfEx.R)
===================================================================
--- pkg/RcppSMC/R/pfLineartBS.R (rev 0)
+++ pkg/RcppSMC/R/pfLineartBS.R 2012-03-18 20:37:55 UTC (rev 3525)
@@ -0,0 +1,50 @@
+
+pfLineartBS<- function(data, particles=1000, plot=FALSE, onlinePlot) {
+
+ # if no data supplied, use default
+ if (missing(data)) data <- getEx1Data()
+
+ if (missing(onlinePlot)) {
+ useOnline <- FALSE
+ onlinePlot <- function() { NULL }
+ } else {
+ useOnline <- TRUE
+ # set up x11
+ x11(width=3,height=3)
+ par(mar=c(3,3,1,1),cex=0.8, pch=19)
+ }
+
+ # more eloquent tests can be added
+ stopifnot(nrow(data) > 0,
+ ncol(data) == 2,
+ colnames(data) == c("x", "y"),
+ class(onlinePlot) == "function")
+
+ res <- .Call("pfLineartBS", as.matrix(data),
+ particles,
+ useOnline,
+ onlinePlot,
+ package="RcppSMC")
+
+ if (plot) {
+ ## plot 5.1 from vignette / paper
+ 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)
+}
+
+pfExOnlinePlot <- function(xm, ym) {
+ plot(xm, ym, xlim=c(-7,0), ylim=c(2,14))
+ Sys.sleep(0.05)
+}
+
Deleted: pkg/RcppSMC/inst/include/pffuncs.h
===================================================================
--- pkg/RcppSMC/inst/include/pffuncs.h 2012-03-18 19:52:08 UTC (rev 3524)
+++ pkg/RcppSMC/inst/include/pffuncs.h 2012-03-18 20:37:55 UTC (rev 3525)
@@ -1,48 +0,0 @@
-// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
-//
-// pffuncs.h: Rcpp wrapper for SMC library -- first example of Johansen (2009)
-//
-// Copyright (C) 2008 Adam Johansen
-// Copyright (C) 2012 Dirk Eddelbuettel and Adam Johansen
-//
-// This file is part of RcppSMC.
-//
-// RcppSMC is free software: you can redistribute it and/or modify it
-// under the terms of the GNU General Public License as published by
-// the Free Software Foundation, either version 2 of the License, or
-// (at your option) any later version.
-//
-// RcppSMC is distributed in the hope that it will be useful, but
-// WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with RcppSMC. If not, see <http://www.gnu.org/licenses/>.
-
-#include "smctc.h"
-
-class cv_state
-{
-public:
- double x_pos, y_pos;
- double x_vel, y_vel;
-};
-
-class cv_obs
-{
-public:
- double x_pos, y_pos;
-};
-
-double logLikelihood(long lTime, const cv_state & X);
-
-smc::particle<cv_state> fInitialise(smc::rng *pRng);
-long fSelect(long lTime, const smc::particle<cv_state> & p, smc::rng *pRng);
-void fMove(long lTime, smc::particle<cv_state> & pFrom, smc::rng *pRng);
-
-extern double nu_x;
-extern double nu_y;
-extern double Delta;
-
-extern std::vector<cv_obs> y;
Copied: pkg/RcppSMC/inst/include/pflineart.h (from rev 3519, pkg/RcppSMC/inst/include/pffuncs.h)
===================================================================
--- pkg/RcppSMC/inst/include/pflineart.h (rev 0)
+++ pkg/RcppSMC/inst/include/pflineart.h 2012-03-18 20:37:55 UTC (rev 3525)
@@ -0,0 +1,48 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
+//
+// pffuncs.h: Rcpp wrapper for SMC library -- first example of Johansen (2009)
+//
+// Copyright (C) 2008 Adam Johansen
+// Copyright (C) 2012 Dirk Eddelbuettel and Adam Johansen
+//
+// This file is part of RcppSMC.
+//
+// RcppSMC is free software: you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 2 of the License, or
+// (at your option) any later version.
+//
+// RcppSMC is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with RcppSMC. If not, see <http://www.gnu.org/licenses/>.
+
+#include "smctc.h"
+
+class cv_state
+{
+public:
+ double x_pos, y_pos;
+ double x_vel, y_vel;
+};
+
+class cv_obs
+{
+public:
+ double x_pos, y_pos;
+};
+
+double logLikelihood(long lTime, const cv_state & X);
+
+smc::particle<cv_state> fInitialise(smc::rng *pRng);
+long fSelect(long lTime, const smc::particle<cv_state> & p, smc::rng *pRng);
+void fMove(long lTime, smc::particle<cv_state> & pFrom, smc::rng *pRng);
+
+extern double nu_x;
+extern double nu_y;
+extern double Delta;
+
+extern std::vector<cv_obs> y;
Deleted: pkg/RcppSMC/man/pfEx.Rd
===================================================================
--- pkg/RcppSMC/man/pfEx.Rd 2012-03-18 19:52:08 UTC (rev 3524)
+++ pkg/RcppSMC/man/pfEx.Rd 2012-03-18 20:37:55 UTC (rev 3525)
@@ -1,60 +0,0 @@
-\name{pfEx}
-\alias{pfEx}
-\alias{pfExOnlinePlot}
-\title{Particle Filter Example}
-\description{
- The \code{pfEx} function provides a simple example for
- \pkg{RcppSMC}. It is based on the first example in \code{SMCTC} and
- the discussion in Section 5.1 of Johansen (2009). A simple 'vehicle
- tracking' problem of 100 observations is solved with 1000 particles.
-
- The \code{pfExOnlinePlot} function provides a simple default
- \sQuote{online} plotting function that is invoked during the
- estimation process.
-}
-\usage{
- pfEx(data, particles=1000, plot=FALSE, onlinePlot)
- pfExOnlinePlot(xm, ym)
-}
-\arguments{
- \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.}
- \item{onlinePlot}{A user-supplied callback function which is called with the
- x and y position vectors during each iteration of the algorithm; see
- pfExOnlinePlot for a simple example.}
- \item{xm}{Vector with x position.}
- \item{ym}{Vector with y position.}
-}
-\value{
- The function returns a \code{data.frame} containing as many rows as in
- the input data, and four columns corresponding to the estimated \eqn{x}{x} and
- \eqn{y}{y} coordinates as well as the estimated velocity in these two directions.
-}
-\details{
- The \code{pfEx} function provides a simple example for
- \pkg{RcppSMC}. It is based on the \code{pf} example in the
- \code{SMCTC} library, and discussed in the Section 5.1 of his
- corresponding paper (Johansen, 2009).
-
- Using the simple \code{pfExOnlinePlot} function illustrates how
- callbacks into R, for example for plotting, can be made during the
- operation of SMC algorithm.
-}
-\references{
- A. M. Johansen. SMCTC: Sequential Monte Carlo in C++.
- Journal of Statistical Software, 30(6):1-41, April
- 2009. \url{http://www.jstatsoft.org/v30/i06/paper}
-}
-\seealso{The SMCTC paper and code at \url{http://www.jstatsoft.org/v30/i06/paper}.}
-\examples{
-\dontrun{
- res <- pfEx(plot=TRUE)
- res <- pfEx(onlinePlot=pfExOnlinePlot)
-}
-}
-\author{Adam M. Johansen for SMCTC, Dirk Eddelbuettel for the RcppSMC wrapper.}
-\keyword{programming}
Copied: pkg/RcppSMC/man/pfLineartBS.Rd (from rev 3519, pkg/RcppSMC/man/pfEx.Rd)
===================================================================
--- pkg/RcppSMC/man/pfLineartBS.Rd (rev 0)
+++ pkg/RcppSMC/man/pfLineartBS.Rd 2012-03-18 20:37:55 UTC (rev 3525)
@@ -0,0 +1,61 @@
+\name{pfLineartBS}
+\alias{pfLineartBS}
+\alias{pfExOnlinePlot}
+\title{Particle Filter Example}
+\description{
+ The \code{pfLineartBS} function provides a simple example for
+ \pkg{RcppSMC}. It is based on the first example in \code{SMCTC} and
+ the discussion in Section 5.1 of Johansen (2009). A simple 'vehicle
+ tracking' problem of 100 observations is solved with 1000 particles.
+
+ The \code{pfExOnlinePlot} function provides a simple default
+ \sQuote{online} plotting function that is invoked during the
+ estimation process.
+}
+\usage{
+ pfLineartBS(data, particles=1000, plot=FALSE, onlinePlot)
+ pfExOnlinePlot(xm, ym)
+}
+\arguments{
+ \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.}
+ \item{onlinePlot}{A user-supplied callback function which is called with the
+ x and y position vectors during each iteration of the algorithm; see
+ pfExOnlinePlot for a simple example.}
+ \item{xm}{Vector with x position.}
+ \item{ym}{Vector with y position.}
+}
+\value{
+ The function returns a \code{data.frame} containing as many rows as in
+ the input data, and four columns corresponding to the estimated \eqn{x}{x} and
+ \eqn{y}{y} coordinates as well as the estimated velocity in these two directions.
+}
+\details{
+ The \code{pfLineartBS} function provides a simple example for
+ \pkg{RcppSMC}. The model is linear with t-distributed innovations.
+ It is based on the \code{pf} example in the
+ \code{SMCTC} library, and discussed in the Section 5.1 of his
+ corresponding paper (Johansen, 2009).
+
+ Using the simple \code{pfExOnlinePlot} function illustrates how
+ callbacks into R, for example for plotting, can be made during the
+ operation of SMC algorithm.
+}
+\references{
+ A. M. Johansen. SMCTC: Sequential Monte Carlo in C++.
+ Journal of Statistical Software, 30(6):1-41, April
+ 2009. \url{http://www.jstatsoft.org/v30/i06/paper}
+}
+\seealso{The SMCTC paper and code at \url{http://www.jstatsoft.org/v30/i06/paper}.}
+\examples{
+\dontrun{
+ res <- pfLineartBS(plot=TRUE)
+ res <- pfLineartBS(onlinePlot=pfExOnlinePlot)
+}
+}
+\author{Adam M. Johansen for SMCTC, Dirk Eddelbuettel for the RcppSMC wrapper.}
+\keyword{programming}
Deleted: pkg/RcppSMC/src/pf.cpp
===================================================================
--- pkg/RcppSMC/src/pf.cpp 2012-03-18 19:52:08 UTC (rev 3524)
+++ pkg/RcppSMC/src/pf.cpp 2012-03-18 20:37:55 UTC (rev 3525)
@@ -1,211 +0,0 @@
-// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
-//
-// pf.cpp: Rcpp wrapper for SMC library -- first example of Johansen (2009)
-//
-// Copyright (C) 2008 Adam Johansen
-// Copyright (C) 2012 Dirk Eddelbuettel and Adam Johansen
-//
-// This file is part of RcppSMC.
-//
-// RcppSMC is free software: you can redistribute it and/or modify it
-// under the terms of the GNU General Public License as published by
-// the Free Software Foundation, either version 2 of the License, or
-// (at your option) any later version.
-//
-// RcppSMC is distributed in the hope that it will be useful, but
-// WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU General Public License
-// along with RcppSMC. If not, see <http://www.gnu.org/licenses/>.
-
-#include <Rcpp.h>
-
-#include "smctc.h"
-#include "pffuncs.h"
-#include "rngR.h"
-
-#include <cstdio>
-#include <cstdlib>
-#include <cstring>
-
-using namespace std;
-
-///The observations
-//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*);
-double integrand_var_x(const cv_state&, void*);
-double integrand_var_y(const cv_state&, void*);
-
-// pf() function callable from R via Rcpp:: essentially the same as main() from pf.cc
-// minor interface change to pass data down as matrix, rather than a filename
-extern "C" SEXP pf(SEXP dataS, SEXP partS, SEXP usefS, SEXP funS) {
-
- long lIterates;
-
- try {
-
- //std::string filename = Rcpp::as<std::string>(fileS);
- unsigned long lNumber = Rcpp::as<unsigned long>(partS);
- bool useF = Rcpp::as<bool>(usefS);
- Rcpp::Function f(funS);
-
- // 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);
- smc::moveset<cv_state> Moveset(fInitialise, fMove, NULL);
-
- Sampler.SetResampleParams(SMC_RESAMPLE_RESIDUAL, 0.5);
- Sampler.SetMoveSet(Moveset);
- Sampler.Initialise();
-
- Rcpp::NumericVector Xm(lIterates), Xv(lIterates), Ym(lIterates), Yv(lIterates);
-
- for(int n=0; n < lIterates; ++n) {
- Sampler.Iterate();
-
- Xm(n) = Sampler.Integrate(integrand_mean_x, NULL);
- Xv(n) = Sampler.Integrate(integrand_var_x, (void*)&Xm(n));
- Ym(n) = Sampler.Integrate(integrand_mean_y, NULL);
- Yv(n) = Sampler.Integrate(integrand_var_y, (void*)&Ym(n));
-
- if (useF) f(Xm, Ym);
- }
-
- return Rcpp::DataFrame::create(Rcpp::Named("Xm") = Xm,
- Rcpp::Named("Xv") = Xv,
- Rcpp::Named("Ym") = Ym,
- Rcpp::Named("Yv") = Yv);
- }
- catch(smc::exception e) {
- Rcpp::Rcout << e; // not cerr, as R doesn't like to mix with i/o
- //exit(e.lCode); // we're just called from R so we should not exit
- }
- return R_NilValue; // to provide a return
-}
-
-#if 0
-long load_data(char const * szName, cv_obs** yp)
-{
- FILE * fObs = fopen(szName,"rt");
- if (!fObs)
- throw SMC_EXCEPTION(SMCX_FILE_NOT_FOUND, "Error: pf assumes that the current directory contains an appropriate data file called data.csv\nThe first line should contain a constant indicating the number of data lines it contains.\nThe remaining lines should contain comma-separated pairs of x,y observations.");
- char szBuffer[1024];
- char* rc = fgets(szBuffer, 1024, fObs);
- if (rc==NULL)
- throw SMC_EXCEPTION(SMCX_FILE_NOT_FOUND, "Error: no data found.");
- long lIterates = strtol(szBuffer, NULL, 10);
-
- *yp = new cv_obs[lIterates];
-
- for(long i = 0; i < lIterates; ++i)
- {
- rc = fgets(szBuffer, 1024, fObs);
- (*yp)[i].x_pos = strtod(strtok(szBuffer, ",\r\n "), NULL);
- (*yp)[i].y_pos = strtod(strtok(szBuffer, ",\r\n "), NULL);
- }
- fclose(fObs);
-
- return lIterates;
-}
-#endif
-
-double integrand_mean_x(const cv_state& s, void *)
-{
- return s.x_pos;
-}
-
-double integrand_var_x(const cv_state& s, void* vmx)
-{
- double* dmx = (double*)vmx;
- double d = (s.x_pos - (*dmx));
- return d*d;
-}
-
-double integrand_mean_y(const cv_state& s, void *)
-{
- return s.y_pos;
-}
-
-double integrand_var_y(const cv_state& s, void* vmy)
-{
- double* dmy = (double*)vmy;
- double d = (s.y_pos - (*dmy));
- return d*d;
-}
-#include <iostream>
-#include <cmath>
-//#include <gsl/gsl_randist.h>
-
-
-using namespace std;
-
-double var_s0 = 4;
-double var_u0 = 1;
-double var_s = 0.02;
-double var_u = 0.001;
-
-double scale_y = 0.1;
-double nu_y = 10.0;
-double Delta = 0.1;
-
-///The function corresponding to the log likelihood at specified time and position (up to normalisation)
-
-/// \param lTime The current time (i.e. the index of the current distribution)
-/// \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));
-
- //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
-
-/// \param pRng A pointer to the random number generator which is to be used
-smc::particle<cv_state> fInitialise(smc::rng *pRng)
-{
- cv_state value;
-
- value.x_pos = pRng->Normal(0,sqrt(var_s0));
- value.y_pos = pRng->Normal(0,sqrt(var_s0));
- value.x_vel = pRng->Normal(0,sqrt(var_u0));
- value.y_vel = pRng->Normal(0,sqrt(var_u0));
-
- return smc::particle<cv_state>(value,logLikelihood(0,value));
-}
-
-///The proposal function.
-
-///\param lTime The sampler iteration.
-///\param pFrom The particle to move.
-///\param pRng A random number generator.
-void fMove(long lTime, smc::particle<cv_state > & pFrom, smc::rng *pRng)
-{
- cv_state * cv_to = pFrom.GetValuePointer();
-
- cv_to->x_pos += cv_to->x_vel * Delta + pRng->Normal(0,sqrt(var_s));
- cv_to->x_vel += pRng->Normal(0,sqrt(var_u));
- cv_to->y_pos += cv_to->y_vel * Delta + pRng->Normal(0,sqrt(var_s));
- cv_to->y_vel += pRng->Normal(0,sqrt(var_u));
- pFrom.AddToLogWeight(logLikelihood(lTime, *cv_to));
-}
Copied: pkg/RcppSMC/src/pflineart.cpp (from rev 3519, pkg/RcppSMC/src/pf.cpp)
===================================================================
--- pkg/RcppSMC/src/pflineart.cpp (rev 0)
+++ pkg/RcppSMC/src/pflineart.cpp 2012-03-18 20:37:55 UTC (rev 3525)
@@ -0,0 +1,211 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
+//
+// pflineart.cpp: Rcpp wrapper for SMC library -- first example of Johansen (2009)
+//
+// Copyright (C) 2008 Adam Johansen
+// Copyright (C) 2012 Dirk Eddelbuettel and Adam Johansen
+//
+// This file is part of RcppSMC.
+//
+// RcppSMC is free software: you can redistribute it and/or modify it
+// under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 2 of the License, or
+// (at your option) any later version.
+//
+// RcppSMC is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with RcppSMC. If not, see <http://www.gnu.org/licenses/>.
+
+#include <Rcpp.h>
+
+#include "smctc.h"
+#include "pflineart.h"
+#include "rngR.h"
+
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+
+using namespace std;
+
+///The observations
+//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*);
+double integrand_var_x(const cv_state&, void*);
+double integrand_var_y(const cv_state&, void*);
+
+// pf() function callable from R via Rcpp:: essentially the same as main() from pf.cc
+// minor interface change to pass data down as matrix, rather than a filename
+extern "C" SEXP pfLineartBS(SEXP dataS, SEXP partS, SEXP usefS, SEXP funS) {
+
+ long lIterates;
+
+ try {
+
+ //std::string filename = Rcpp::as<std::string>(fileS);
+ unsigned long lNumber = Rcpp::as<unsigned long>(partS);
+ bool useF = Rcpp::as<bool>(usefS);
+ Rcpp::Function f(funS);
+
+ // 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);
+ smc::moveset<cv_state> Moveset(fInitialise, fMove, NULL);
+
+ Sampler.SetResampleParams(SMC_RESAMPLE_RESIDUAL, 0.5);
+ Sampler.SetMoveSet(Moveset);
+ Sampler.Initialise();
+
+ Rcpp::NumericVector Xm(lIterates), Xv(lIterates), Ym(lIterates), Yv(lIterates);
+
+ for(int n=0; n < lIterates; ++n) {
+ Sampler.Iterate();
+
+ Xm(n) = Sampler.Integrate(integrand_mean_x, NULL);
+ Xv(n) = Sampler.Integrate(integrand_var_x, (void*)&Xm(n));
+ Ym(n) = Sampler.Integrate(integrand_mean_y, NULL);
+ Yv(n) = Sampler.Integrate(integrand_var_y, (void*)&Ym(n));
+
+ if (useF) f(Xm, Ym);
+ }
+
+ return Rcpp::DataFrame::create(Rcpp::Named("Xm") = Xm,
+ Rcpp::Named("Xv") = Xv,
+ Rcpp::Named("Ym") = Ym,
+ Rcpp::Named("Yv") = Yv);
+ }
+ catch(smc::exception e) {
+ Rcpp::Rcout << e; // not cerr, as R doesn't like to mix with i/o
+ //exit(e.lCode); // we're just called from R so we should not exit
+ }
+ return R_NilValue; // to provide a return
+}
+
+#if 0
+long load_data(char const * szName, cv_obs** yp)
+{
+ FILE * fObs = fopen(szName,"rt");
+ if (!fObs)
+ throw SMC_EXCEPTION(SMCX_FILE_NOT_FOUND, "Error: pf assumes that the current directory contains an appropriate data file called data.csv\nThe first line should contain a constant indicating the number of data lines it contains.\nThe remaining lines should contain comma-separated pairs of x,y observations.");
+ char szBuffer[1024];
+ char* rc = fgets(szBuffer, 1024, fObs);
+ if (rc==NULL)
+ throw SMC_EXCEPTION(SMCX_FILE_NOT_FOUND, "Error: no data found.");
+ long lIterates = strtol(szBuffer, NULL, 10);
+
+ *yp = new cv_obs[lIterates];
+
+ for(long i = 0; i < lIterates; ++i)
+ {
+ rc = fgets(szBuffer, 1024, fObs);
+ (*yp)[i].x_pos = strtod(strtok(szBuffer, ",\r\n "), NULL);
+ (*yp)[i].y_pos = strtod(strtok(szBuffer, ",\r\n "), NULL);
+ }
+ fclose(fObs);
+
+ return lIterates;
+}
+#endif
+
+double integrand_mean_x(const cv_state& s, void *)
+{
+ return s.x_pos;
+}
+
+double integrand_var_x(const cv_state& s, void* vmx)
+{
+ double* dmx = (double*)vmx;
+ double d = (s.x_pos - (*dmx));
+ return d*d;
+}
+
+double integrand_mean_y(const cv_state& s, void *)
+{
+ return s.y_pos;
+}
+
+double integrand_var_y(const cv_state& s, void* vmy)
+{
+ double* dmy = (double*)vmy;
+ double d = (s.y_pos - (*dmy));
+ return d*d;
+}
+#include <iostream>
+#include <cmath>
+//#include <gsl/gsl_randist.h>
+
+
+using namespace std;
+
+double var_s0 = 4;
+double var_u0 = 1;
+double var_s = 0.02;
+double var_u = 0.001;
+
+double scale_y = 0.1;
+double nu_y = 10.0;
+double Delta = 0.1;
+
+///The function corresponding to the log likelihood at specified time and position (up to normalisation)
+
+/// \param lTime The current time (i.e. the index of the current distribution)
+/// \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));
+
+ //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
+
+/// \param pRng A pointer to the random number generator which is to be used
+smc::particle<cv_state> fInitialise(smc::rng *pRng)
+{
+ cv_state value;
+
+ value.x_pos = pRng->Normal(0,sqrt(var_s0));
+ value.y_pos = pRng->Normal(0,sqrt(var_s0));
+ value.x_vel = pRng->Normal(0,sqrt(var_u0));
+ value.y_vel = pRng->Normal(0,sqrt(var_u0));
+
+ return smc::particle<cv_state>(value,logLikelihood(0,value));
+}
+
+///The proposal function.
+
+///\param lTime The sampler iteration.
+///\param pFrom The particle to move.
+///\param pRng A random number generator.
+void fMove(long lTime, smc::particle<cv_state > & pFrom, smc::rng *pRng)
+{
+ cv_state * cv_to = pFrom.GetValuePointer();
+
+ cv_to->x_pos += cv_to->x_vel * Delta + pRng->Normal(0,sqrt(var_s));
+ cv_to->x_vel += pRng->Normal(0,sqrt(var_u));
+ cv_to->y_pos += cv_to->y_vel * Delta + pRng->Normal(0,sqrt(var_s));
+ cv_to->y_vel += pRng->Normal(0,sqrt(var_u));
+ pFrom.AddToLogWeight(logLikelihood(lTime, *cv_to));
+}
More information about the Rcpp-commits
mailing list