[Rcpp-commits] r3511 - in pkg/RcppSMC: . R inst/include man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 16 12:09:11 CET 2012
Author: adamj
Date: 2012-03-16 12:09:11 +0100 (Fri, 16 Mar 2012)
New Revision: 3511
Added:
pkg/RcppSMC/R/pfNonlinBS.R
pkg/RcppSMC/inst/include/pfnonlinbs.h
pkg/RcppSMC/man/pfNonlinBS.Rd
pkg/RcppSMC/src/pfnonlinbs.cpp
Modified:
pkg/RcppSMC/ChangeLog
pkg/RcppSMC/NAMESPACE
Log:
Added simple PF example from Gordon, Salmond and Smith (1993) to RcppSMC.
Modified: pkg/RcppSMC/ChangeLog
===================================================================
--- pkg/RcppSMC/ChangeLog 2012-03-15 21:03:58 UTC (rev 3510)
+++ pkg/RcppSMC/ChangeLog 2012-03-16 11:09:11 UTC (rev 3511)
@@ -1,3 +1,11 @@
+2012-03-16 Adam Johansen <a.m.johansen at warwick.ac.uk>
+ * NAMESPACE: Added pfNonlinBS.
+ * src/pfnonlinbS.cpp: Bootstrap particle nonlinear particle filter
+ example.
+ * inst/include/pfnonlinbs.h: Header for the same.
+ * R/pfNonlinBS.R: Minimal wrapper and output for the same.
+ * man/pfNonlinBS.Rd: Minimal documentation for the same.
+
2012-03-15 Dirk Eddelbuettel <edd at debian.org>
* src/pf.cpp: Another new/delete cleanup
Modified: pkg/RcppSMC/NAMESPACE
===================================================================
--- pkg/RcppSMC/NAMESPACE 2012-03-15 21:03:58 UTC (rev 3510)
+++ pkg/RcppSMC/NAMESPACE 2012-03-16 11:09:11 UTC (rev 3511)
@@ -3,4 +3,5 @@
export(blockpfGaussianOpt,
pfEx,
+ pfNonlinBS,
rareEventsEx)
Added: pkg/RcppSMC/R/pfNonlinBS.R
===================================================================
--- pkg/RcppSMC/R/pfNonlinBS.R (rev 0)
+++ pkg/RcppSMC/R/pfNonlinBS.R 2012-03-16 11:09:11 UTC (rev 3511)
@@ -0,0 +1,20 @@
+pfNonlinBS <- function(data=c(), particles=500, plot=FALSE) {
+ if (length(data == 0)) {
+ #Include some error handling here
+ return;
+ }
+ res <- .Call("pfNonlinBS", data, particles, package="RcppSMC")
+
+ time <- 1:length(data);
+ if (plot) {
+ with(res, plot(time,mean,type='l',xlab='time',ylab='estimate'))
+ with(res,polygon(c(time,51-time),c(mean-2*sd,(mean+2*sd)[seq(length(data),1,-1)]),col='palegreen1',border=NA))
+ with(res,polygon(c(time,51-time),c(mean-1*sd,(mean+1*sd)[seq(length(data),1,-1)]),col='palegreen3',border=NA))
+ with(res,lines(time,mean, lwd=2, col='darkblue'))
+ }
+
+ invisible(res)
+}
+
+
+
Added: pkg/RcppSMC/inst/include/pfnonlinbs.h
===================================================================
--- pkg/RcppSMC/inst/include/pfnonlinbs.h (rev 0)
+++ pkg/RcppSMC/inst/include/pfnonlinbs.h 2012-03-16 11:09:11 UTC (rev 3511)
@@ -0,0 +1,48 @@
+//
+// pf1dnonlinbs.h
+//
+// The declarations and externals for an implementation of the bootstrap
+// particle filter of "Novel approaches to nonlinear non-Gaussian
+// Bayesian state estimation", Gordon Salmond and Smith,
+// IEE PROCEEDINGS-F 140(2):107-113, 1993
+//
+// Copyright Adam Johansen, 2012
+//
+// 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 RInside. If not, see <http://www.gnu.org/licenses/>.
+// from examples/pf/pfexample.cc; pffuncs.cc and pffuncs.hh also used
+
+// RcppSMC builds on top of, and wraps, SMCTC which is
+//
+// Copyright Adam Johansen, 2008-2012.
+//
+// and released under GPL-3, see the copyright headers in inst/include/
+
+
+#include "smctc.h"
+
+namespace nonlinbs {
+ double logLikelihood(long lTime, const double & X);
+
+ smc::particle<double> fInitialise(smc::rng *pRng);
+ long fSelect(long lTime, const smc::particle<double> & p,
+ smc::rng *pRng);
+ void fMove(long lTime, smc::particle<double> & pFrom,
+ smc::rng *pRng);
+
+ double integrand_mean_x(const double&, void*);
+ double integrand_var_x(const double&, void*);
+}
Added: pkg/RcppSMC/man/pfNonlinBS.Rd
===================================================================
--- pkg/RcppSMC/man/pfNonlinBS.Rd (rev 0)
+++ pkg/RcppSMC/man/pfNonlinBS.Rd 2012-03-16 11:09:11 UTC (rev 3511)
@@ -0,0 +1,46 @@
+\name{pfNonlinBS}
+\alias{pfNonlinBS}
+\title{Block Sampling Particle Filter (Linear Gaussian Model; Optimal Proposal)}
+\description{
+ The \code{blockpfGaussianOpt} function provides a simple example for
+ \pkg{RcppSMC}. It is a simple "bootstrap" particle filter which employs
+ multinomial resampling after each iteration applied to the common nonlinear
+ state space model following Gordon, Salmond and Smith (1993).
+}
+\usage{
+ pfNonlinBS(data, particles=500, plot=FALSE)
+}
+\arguments{
+ \item{data}{A vector variable containing the sequence of observations.}
+ \item{particles}{An integer specifying the number of particles.}
+ \item{plot}{A boolean variable describing whether a plot should
+ illustrate the (posterior mean) estimated path along with one and two
+ standard deviation intervals.}
+}
+\value{
+ The function returns two vectors, the first containing the posterior
+ filtering means; the second the posterior filtering standard deviations.
+}
+\details{
+ The \code{pfNonlinbs} function provides a simple example for
+ \pkg{RcppSMC}. It is based on a simple nonlinear state space model in
+ which the state evolution and observation equations are:
+ x(n) = 0.5 x(n-1) + 25 x(n-1) / (1+x(n-1)^2) + 8 cos(1.2(n-1))+ e(n) and
+ y(n) = x(n)^2 / 20 + f(n)
+ where e(n) and f(n) are mutually-independent normal random
+ variables of variances 10.0 and 1.0, respectively. A boostrap proposal
+ (i.e. sampling from the state equation) is used, together with multinomial
+ resampling after each iteration.
+}
+\references{
+ N. J. Gordon, S. J. Salmond, and A. F. M. Smith. Novel approach to
+ nonlinear/non-Gaussian Bayesian state estimation. IEE Proceedings-F,
+ 140(2):107-113, April 1993.
+}
+\examples{
+\dontrun{
+ res <- pfNonlinBS(data,particles=500,plot=TRUE)
+}
+}
+\author{Adam M. Johansen}
+\keyword{programming}
Added: pkg/RcppSMC/src/pfnonlinbs.cpp
===================================================================
--- pkg/RcppSMC/src/pfnonlinbs.cpp (rev 0)
+++ pkg/RcppSMC/src/pfnonlinbs.cpp 2012-03-16 11:09:11 UTC (rev 3511)
@@ -0,0 +1,135 @@
+//
+// pf1dnonlinbs.cpp
+//
+// The principle source file for an implementation of the bootstrap
+// particle filter of "Novel approaches to nonlinear non-Gaussian
+// Bayesian state estimation", Gordon Salmond and Smith,
+// IEE PROCEEDINGS-F 140(2):107-113, 1993
+//
+// Copyright Adam Johansen, 2012
+//
+// 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 RInside. If not, see <http://www.gnu.org/licenses/>.
+// from examples/pf/pfexample.cc; pffuncs.cc and pffuncs.hh also used
+
+// RcppSMC builds on top of, and wraps, SMCTC which is
+//
+// Copyright Adam Johansen, 2008-2012.
+//
+// and released under GPL-3, see the copyright headers in inst/include/
+
+#include "smctc.h"
+#include "pfnonlinbs.h"
+
+#include <cstdlib>
+#include <cmath>
+
+namespace nonlinbs {
+const double std_x0 = 2;
+const double var_x0 = std_x0 * std_x0;
+const double std_x = sqrt(10.0);
+const double var_x = std_x * std_x;
+const double var_y = 1.0;
+
+const double scale_y = 1.0 / 20.0;
+
+///The observations
+Rcpp::NumericVector y;
+}
+
+using namespace std;
+using namespace nonlinbs;
+
+extern "C" SEXP pfNonlinBS(SEXP dataS, SEXP partS)
+{
+ long lNumber = Rcpp::as<long>(partS);
+
+ y = Rcpp::NumericVector(dataS);
+ long lIterates = y.size();
+
+ //Initialise and run the sampler
+ smc::sampler<double> Sampler(lNumber, SMC_HISTORY_NONE);
+ smc::moveset<double> Moveset(fInitialise, fMove, NULL);
+
+ Sampler.SetResampleParams(SMC_RESAMPLE_MULTINOMIAL, 1.01 * lNumber);
+ Sampler.SetMoveSet(Moveset);
+ Sampler.Initialise();
+
+ Rcpp::NumericVector resMean = Rcpp::NumericVector(lIterates);
+ Rcpp::NumericVector resSD = Rcpp::NumericVector(lIterates);
+ for(int n=0 ; n < lIterates ; ++n) {
+ if(n > 0)
+ Sampler.Iterate();
+
+ resMean(n) = Sampler.Integrate(integrand_mean_x,NULL);
+ resSD(n) = sqrt(Sampler.Integrate(integrand_var_x, (void*)&resSD(n)));
+ }
+
+ return Rcpp::List::create(Rcpp::_["mean"] = resMean,
+ Rcpp::_["sd"] = resSD);
+}
+
+namespace nonlinbs {
+///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 double & x)
+{
+ return -0.5 * pow(y[int(lTime)] - x*x*scale_y,2) / var_y;
+}
+
+///A function to initialise particles
+
+/// \param pRng A pointer to the random number generator which is to be used
+smc::particle<double> fInitialise(smc::rng *pRng)
+{
+ double x;
+
+ x = pRng->Normal(0,std_x0);
+
+ return smc::particle<double>(x,logLikelihood(0,x));
+}
+
+///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<double> & pFrom, smc::rng *pRng)
+{
+ double *to = pFrom.GetValuePointer();
+
+ double x = 0.5 * (*to) + 25.0*(*to) / (1.0 + (*to) * (*to)) + 8.0 * cos(1.2 * ( lTime)) + pRng->Normal(0.0,std_x);
+
+ *to = x;
+
+ pFrom.AddToLogWeight(logLikelihood(lTime, *to));
+}
+
+
+double integrand_mean_x(const double& x, void *)
+{
+ return x;
+}
+
+double integrand_var_x(const double& x, void* vmx)
+{
+ double* dmx = (double*)vmx;
+ double d = (x - (*dmx));
+ return d*d;
+}
+}
More information about the Rcpp-commits
mailing list