[Rcpp-commits] r3502 - in pkg/RcppSMC: . R inst/include man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Mar 13 18:49:37 CET 2012
Author: adamj
Date: 2012-03-13 18:49:37 +0100 (Tue, 13 Mar 2012)
New Revision: 3502
Added:
pkg/RcppSMC/R/bspfGaussianOptimal.R
pkg/RcppSMC/inst/include/blockpfgaussianopt.hh
pkg/RcppSMC/man/blockpfGaussianOpt.Rd
pkg/RcppSMC/src/blockpfgaussianopt.cc
pkg/RcppSMC/src/blockpfgaussianoptfuncs.cc
Modified:
pkg/RcppSMC/ChangeLog
pkg/RcppSMC/DESCRIPTION
pkg/RcppSMC/NAMESPACE
Log:
Adding a block-sampling particle filter example to RcppSMC. This toy example
allows data to be passed from R to SMCTC and passes results back again.
Modified: pkg/RcppSMC/ChangeLog
===================================================================
--- pkg/RcppSMC/ChangeLog 2012-03-10 20:20:20 UTC (rev 3501)
+++ pkg/RcppSMC/ChangeLog 2012-03-13 17:49:37 UTC (rev 3502)
@@ -1,3 +1,17 @@
+2012-03-13 Adam Johansen <a.m.johansen at warwick.ac.uk>
+
+ * DESCRIPTION: Version 0.0.3
+
+ * NAMESPACE: Added blockpfGaussianOpt.
+ * inst/include/blockpfguassianopt.hh: Declarations for univariate
+ block-sampling particle filter implementation.
+ * src/blockpfguassianopt.cc: Main file for univariate block-sampling
+ particle filter implementation.
+ * src/blockpfguassianoptfuncs.cc: Internal functions for same.
+ * R/blockpfguassianopt.R: Minimal wrapper and output plotting for
+ same.
+ * doc/blockpfguassianopt.Rd: Added minimal documentation.
+
2012-01-21 Dirk Eddelbuettel <edd at debian.org>
* DESCRIPTION: Version 0.0.2
Modified: pkg/RcppSMC/DESCRIPTION
===================================================================
--- pkg/RcppSMC/DESCRIPTION 2012-03-10 20:20:20 UTC (rev 3501)
+++ pkg/RcppSMC/DESCRIPTION 2012-03-13 17:49:37 UTC (rev 3502)
@@ -1,17 +1,17 @@
Package: RcppSMC
Type: Package
-Title: Rcpp bindings for Sequential Monte Carlo
-Version: 0.0.2
+Title: Rcpp bindings for Sequential Monte Carlo
+Version: 0.0.3
Date: $Date$
-Author: Dirk Eddelbuettel
+Author: Dirk Eddelbuettel and Adam M. Johansen
Maintainer: Dirk Eddelbuettel <edd at debian.org>
Description: This package provides R with access to the Sequential Monte
Carlo Template Classes by Johansen (J Stat Soft, 2009, v30, i6).
.
- At present, this is a proof of concept. SMCTC uses GSL RNGs which we can
- easily swap for R RNGs. Further integration pending.
+ At present, this is a proof of concept.
+ Further integration pending.
License: GPL (>= 2)
LazyLoad: yes
Depends: Rcpp (>= 0.9.0), methods
LinkingTo: Rcpp
-
+Packaged: 2012-03-13 17:50:23 UTC; adamj
Modified: pkg/RcppSMC/NAMESPACE
===================================================================
--- pkg/RcppSMC/NAMESPACE 2012-03-10 20:20:20 UTC (rev 3501)
+++ pkg/RcppSMC/NAMESPACE 2012-03-13 17:49:37 UTC (rev 3502)
@@ -1,5 +1,6 @@
useDynLib(RcppSMC)
-export(pfEx,
+export(blockpfGaussianOpt,
+ pfEx,
rareEventsEx)
Added: pkg/RcppSMC/R/bspfGaussianOptimal.R
===================================================================
--- pkg/RcppSMC/R/bspfGaussianOptimal.R (rev 0)
+++ pkg/RcppSMC/R/bspfGaussianOptimal.R 2012-03-13 17:49:37 UTC (rev 3502)
@@ -0,0 +1,32 @@
+
+blockpfGaussianOpt <- function(data=c(), particles=1000, lag=5, plot=FALSE) {
+
+ if (length(data == 0)) {
+ #Include some error handling here
+ return;
+ }
+ res <- .Call("blockpfGaussianOpt", data, particles, lag, package="RcppSMC")
+
+ if (plot) {
+ time = 1:length(data);
+ mvect = t(res$weight) %*% res$values / sum(res$weight);
+ sqvect = t(res$weight) %*% res$values^2 / sum(res$weight);
+ sdvect = sqrt(sqvect - mvect^2);
+
+ plot(time, mvect, col='dark red', 'l', lty = 1, xlab = 'Iteration',
+ ylab='State', main='Mean and 1, 2 standard deviation credible
+ intervals with observations', xlim = c(0,length(data)), ylim=c(min(mvect - 2.1
+ * (sdvect)), max(mvect+2.1*sdvect))
+ )
+ lines(time, mvect + sdvect, lty=3, col='dark blue')
+ lines(time, mvect - sdvect, lty=3, col='dark blue')
+ lines(time, mvect + 2 * sdvect, lty=2, col='dark blue')
+ lines(time, mvect - 2 * sdvect, lty=2, col='dark blue')
+ points(time, data, col = 'dark green')
+ }
+
+ invisible(res)
+}
+
+
+
Added: pkg/RcppSMC/inst/include/blockpfgaussianopt.hh
===================================================================
--- pkg/RcppSMC/inst/include/blockpfgaussianopt.hh (rev 0)
+++ pkg/RcppSMC/inst/include/blockpfgaussianopt.hh 2012-03-13 17:49:37 UTC (rev 3502)
@@ -0,0 +1,15 @@
+#include "smctc.hh"
+#include <vector>
+
+using namespace std;
+
+smc::particle<vector<double> > fInitialiseBSPFG(smc::rng *pRng);
+//long fSelectBSPFG(long lTime, const smc::particle<vector<double> > & p, smc::rng *pRng);
+void fMoveBSPFG(long lTime, smc::particle<vector<double> > & pFrom, smc::rng *pRng);
+
+namespace BSPFG {
+extern Rcpp::NumericVector y;
+}
+using BSPFG::y;
+
+extern long lLag;
Added: pkg/RcppSMC/man/blockpfGaussianOpt.Rd
===================================================================
--- pkg/RcppSMC/man/blockpfGaussianOpt.Rd (rev 0)
+++ pkg/RcppSMC/man/blockpfGaussianOpt.Rd 2012-03-13 17:49:37 UTC (rev 3502)
@@ -0,0 +1,46 @@
+\name{blockpfGaussianOpt}
+\alias{blockpfGaussianOpt}
+\title{Block Sampling Particle Filter (Linear Gaussian Model; Optimal Proposal)}
+\description{
+ The \code{blockpfGaussianOpt} function provides a simple example for
+ \pkg{RcppSMC}. It is based on a block sampling particle filter for a linear
+ Gaussian model. This is intended only to illustrate the potential of block
+ sampling; one would not ordinarily use a particle filter for a model in
+ which analytic solutions are available. The 'optimal' block sampler in the
+ sense of (Doucet, Briers and Senecal, 2006) can be implemented in this case.
+}
+\usage{
+ blockpfGaussianOpt(data, particles=1000, lag=5, plot=FALSE)
+}
+\arguments{
+ \item{data}{A vector variable containing the sequence of observations.}
+ \item{particles}{An integer specifying the number of particles.}
+ \item{lag}{An integer specifying the length of block to use.}
+ \item{plot}{A boolean variable describing whether plot should
+ illustrate the estimated path along with the uncertainty.}
+}
+\value{
+ The function returns a matrix containing the final sample paths and a vector
+ containing their weights.
+}
+\details{
+ The \code{pfEx} function provides a simple example for
+ \pkg{RcppSMC}. It is based on a simple linear Gaussian state space model in
+ which the state evolution and observation equations are:
+ x(n) = x(n-1) + e(n) and
+ y(n) = x(n) + f(n)
+ where e(n) and f(n) are mutually-independent standard normal random variables. The 'optimal'
+ block-sampling proposal described by (Doucet et al., 2006) is employed.
+}
+\references{
+ A. Doucet, M. Briers, and S. Senecal. Efficient Block Sampling Strategies
+ for sequential Monte Carlo methods. Journal of Computational and Graphical
+ Statistics, 15(3):693-711, 2006.
+}
+\examples{
+\dontrun{
+ res <- blockpfGaussianOpt(data,lag=5,plot=TRUE)
+}
+}
+\author{Adam M. Johansen}
+\keyword{programming}
Added: pkg/RcppSMC/src/blockpfgaussianopt.cc
===================================================================
--- pkg/RcppSMC/src/blockpfgaussianopt.cc (rev 0)
+++ pkg/RcppSMC/src/blockpfgaussianopt.cc 2012-03-13 17:49:37 UTC (rev 3502)
@@ -0,0 +1,59 @@
+#include "smctc.hh"
+#include "blockpfgaussianopt.hh"
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <vector>
+#include <string>
+
+using namespace std;
+
+///The observations
+//double *y;
+namespace BSPFG {
+Rcpp::NumericVector y;
+};
+
+long lLag = 1;
+long load_data(char const * szName, double** y);
+
+extern "C" SEXP blockpfGaussianOpt(SEXP dataS, SEXP partS, SEXP lagS)
+//int main(int argc, char** argv)
+{
+ //Load observations
+ long lIterates;
+ long lNumber = Rcpp::as<long>(partS);
+ lLag = Rcpp::as<long>(lagS);
+
+ y = Rcpp::NumericVector(dataS);
+ lIterates = y.size();
+
+ //Initialise and run the sampler
+ smc::sampler<vector<double> > Sampler(lNumber, SMC_HISTORY_NONE);
+ smc::moveset<vector<double> > Moveset(fInitialiseBSPFG, fMoveBSPFG, NULL);
+
+ Sampler.SetResampleParams(SMC_RESAMPLE_SYSTEMATIC, 0.5);
+ Sampler.SetMoveSet(Moveset);
+
+ Sampler.Initialise();
+ Sampler.IterateUntil(lIterates - 1);
+
+ //Generate results
+ Rcpp::NumericMatrix resValues = Rcpp::NumericMatrix(lNumber,lIterates);
+ Rcpp::NumericVector resWeights = Rcpp::NumericVector(lNumber);
+ for(int i = 0; i < lNumber; ++i)
+ {
+ vector<double> pValue = Sampler.GetParticleValue(i);
+ for(int j = 0; j < lIterates; ++j) {
+ resValues(i,j) = pValue.at(j);
+ }
+ resWeights(i) = Sampler.GetParticleWeight(i);
+ }
+
+ return Rcpp::List::create(Rcpp::_["weight"] = resWeights, Rcpp::_["values"] = resValues);
+}
+
Added: pkg/RcppSMC/src/blockpfgaussianoptfuncs.cc
===================================================================
--- pkg/RcppSMC/src/blockpfgaussianoptfuncs.cc (rev 0)
+++ pkg/RcppSMC/src/blockpfgaussianoptfuncs.cc 2012-03-13 17:49:37 UTC (rev 3502)
@@ -0,0 +1,73 @@
+#include <iostream>
+#include <cmath>
+//#include <gsl/gsl_randist.h>
+
+#include "rngR.h"
+#include "smctc.hh"
+#include "blockpfgaussianopt.hh"
+
+using namespace std;
+using BSPFG::y;
+
+/// \param pRng A pointer to the random number generator which is to be used
+smc::particle<vector<double> > fInitialiseBSPFG(smc::rng *pRng)
+{
+ vector<double> value;
+
+ value.push_back(pRng->Normal(0.5 * y[0],1.0/sqrt(2.0)));
+
+ return smc::particle<vector<double> >(value,1.0);
+}
+
+///The proposal function.
+
+///\param lTime The sampler iteration.
+///\param pFrom The particle to move.
+///\param pRng A random number generator.
+void fMoveBSPFG(long lTime, smc::particle<vector<double> > & pFrom, smc::rng *pRng)
+{
+ std::vector<double> * cv_to = pFrom.GetValuePointer();
+
+ if(lTime == 1) {
+ cv_to->push_back((cv_to->at(lTime-1) + y[int(lTime)])/2.0 + pRng->Normal(0.0,1.0/sqrt(2.0)));
+
+ pFrom.AddToLogWeight(-0.25*(y[int(lTime)] - cv_to->at(lTime-1))*(y[int(lTime)]-cv_to->at(lTime-1)));
+
+ return;
+ }
+
+ long lag = min(lTime,lLag);
+
+ //These structures should really be made static
+ double* mu = new double[lag+1];
+ double* sigma = new double[lag+1];
+ double* sigmah = new double[lag+1];
+ double* mub = new double[lag+1];
+
+ // Forward filtering
+ mu[0] = cv_to->at(lTime-lag);
+ sigma[0] = 0;
+ for(int i = 1; i <= lag; ++i)
+ {
+ sigmah[i] = sigma[i-1] + 1;
+
+ mu[i] = (sigmah[i] * y[int(lTime-lag+i)] + mu[i-1]) / (sigmah[i] + 1);
+ sigma[i] = sigmah[i] / (sigmah[i] + 1);
+ }
+ // Backward smoothing
+ mub[lag] = mu[lag];
+ cv_to->push_back(pRng->Normal(mub[lag],sqrt(sigma[lag])));
+ for(int i = lag-1; i; --i)
+ {
+ mub[i] = (sigma[i]*cv_to->at(lTime-lag+i+1) + mu[i]) / (sigma[i]+1);
+ cv_to->at(lTime-lag+i) = pRng->Normal(mub[i],sqrt(sigma[lag]/(sigma[lag] + 1)));
+ }
+
+ // Importance weighting
+ pFrom.AddToLogWeight(-0.5 * pow(y[int(lTime)] - mu[lag-1],2.0) / (sigmah[lag]+1) );
+
+ delete [] mu;
+ delete [] sigma;
+ delete [] sigmah;
+ delete [] mub;
+}
More information about the Rcpp-commits
mailing list