[Rcpp-commits] r3515 - in pkg/RcppSMC: . R R/rareEvents inst/include inst/include/rareEvents man man/rareEventsEx src src/rareEvents
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Mar 17 10:47:08 CET 2012
Author: adamj
Date: 2012-03-17 10:47:08 +0100 (Sat, 17 Mar 2012)
New Revision: 3515
Added:
pkg/RcppSMC/R/rareEvents/
pkg/RcppSMC/R/rareEvents/rareEventsEx.R
pkg/RcppSMC/inst/include/rareEvents/
pkg/RcppSMC/inst/include/rareEvents/markovchains/
pkg/RcppSMC/inst/include/rareEvents/simfunctions.h
pkg/RcppSMC/man/rareEventsEx/
pkg/RcppSMC/man/rareEventsEx/rareEventsEx.Rd
pkg/RcppSMC/src/rareEvents/
pkg/RcppSMC/src/rareEvents/rareEvents.cpp
pkg/RcppSMC/src/rareEvents/simfunctions.cpp
Removed:
pkg/RcppSMC/R/rareEventsEx.R
pkg/RcppSMC/inst/include/markovchains/
pkg/RcppSMC/inst/include/simfunctions.h
pkg/RcppSMC/man/rareEventsEx.Rd
pkg/RcppSMC/src/rareEvents.cpp
pkg/RcppSMC/src/simfunctions.cpp
Modified:
pkg/RcppSMC/.Rbuildignore
pkg/RcppSMC/ChangeLog
pkg/RcppSMC/NAMESPACE
pkg/RcppSMC/TODO
Log:
Moving rareEvents example.
Modified: pkg/RcppSMC/.Rbuildignore
===================================================================
--- pkg/RcppSMC/.Rbuildignore 2012-03-17 02:27:01 UTC (rev 3514)
+++ pkg/RcppSMC/.Rbuildignore 2012-03-17 09:47:08 UTC (rev 3515)
@@ -1 +1,4 @@
src/compareGSLtoR
+src/rareEvents
+inst/include/rareEvents
+man/rareEvents
Modified: pkg/RcppSMC/ChangeLog
===================================================================
--- pkg/RcppSMC/ChangeLog 2012-03-17 02:27:01 UTC (rev 3514)
+++ pkg/RcppSMC/ChangeLog 2012-03-17 09:47:08 UTC (rev 3515)
@@ -1,3 +1,13 @@
+2012-03-17 Adam Johansen <a.m.johansen at warwick.ac.uk>
+
+ * NAMESPACE removed rareEvents
+ * src/rareEvents.cpp moved to src/rareEvents.
+ * src/simfunctions.cpp idem.
+ * inst/include/simfunctions.h moved to inst/include/rareEvents.
+ * inst/include/markovchains idem
+ * man/rareEventsEx.Rd moved to man/rareEvents.
+ * R/rareEventsEx.R moved to R/rareEvents.
+
2012-03-16 Dirk Eddelbuettel <edd at dexter>
* src/pf.cpp:
Modified: pkg/RcppSMC/NAMESPACE
===================================================================
--- pkg/RcppSMC/NAMESPACE 2012-03-17 02:27:01 UTC (rev 3514)
+++ pkg/RcppSMC/NAMESPACE 2012-03-17 09:47:08 UTC (rev 3515)
@@ -4,5 +4,5 @@
export(blockpfGaussianOpt,
pfEx,
pfExOnlinePlot,
- pfNonlinBS,
- rareEventsEx)
+ pfNonlinBS)
+
Copied: pkg/RcppSMC/R/rareEvents/rareEventsEx.R (from rev 3514, pkg/RcppSMC/R/rareEventsEx.R)
===================================================================
--- pkg/RcppSMC/R/rareEvents/rareEventsEx.R (rev 0)
+++ pkg/RcppSMC/R/rareEvents/rareEventsEx.R 2012-03-17 09:47:08 UTC (rev 3515)
@@ -0,0 +1,9 @@
+## example 2 of Johansen (2009)
+
+rareEventsEx <- function(number=100, iterations=10, threshold=5.0, schedule=30.0) {
+
+ res <- .Call("rareEvents", number, iterations, threshold, schedule, package="RcppSMC")
+
+ invisible(res)
+}
+
Deleted: pkg/RcppSMC/R/rareEventsEx.R
===================================================================
--- pkg/RcppSMC/R/rareEventsEx.R 2012-03-17 02:27:01 UTC (rev 3514)
+++ pkg/RcppSMC/R/rareEventsEx.R 2012-03-17 09:47:08 UTC (rev 3515)
@@ -1,9 +0,0 @@
-## example 2 of Johansen (2009)
-
-rareEventsEx <- function(number=100, iterations=10, threshold=5.0, schedule=30.0) {
-
- res <- .Call("rareEvents", number, iterations, threshold, schedule, package="RcppSMC")
-
- invisible(res)
-}
-
Modified: pkg/RcppSMC/TODO
===================================================================
--- pkg/RcppSMC/TODO 2012-03-17 02:27:01 UTC (rev 3514)
+++ pkg/RcppSMC/TODO 2012-03-17 09:47:08 UTC (rev 3515)
@@ -1,8 +1,6 @@
* Add more examples
- * Maybe remove example 2 from the JSS paper
-
* Maybe rename example 1 from the JSS paper -- what would be a good
name instead of 'pfEx'? Idem for the C++ code.
Copied: pkg/RcppSMC/inst/include/rareEvents/simfunctions.h (from rev 3514, pkg/RcppSMC/inst/include/simfunctions.h)
===================================================================
--- pkg/RcppSMC/inst/include/rareEvents/simfunctions.h (rev 0)
+++ pkg/RcppSMC/inst/include/rareEvents/simfunctions.h 2012-03-17 09:47:08 UTC (rev 3515)
@@ -0,0 +1,32 @@
+#include "smctc.h"
+
+#include "markovchains/markovchain.h"
+
+extern long lIterates;
+extern long lNumber;
+extern long lChainLength;
+extern double dSchedule;
+extern double dThreshold;
+
+double logDensity(long lTime, const mChain<double> & X);
+
+smc::particle<mChain<double> > fInitialiseMC(smc::rng *pRng);
+long fSelect(long lTime, const smc::particle<mChain<double> > & p, smc::rng *pRng);
+void fMove1(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng);
+void fMove2(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng);
+int fMCMC(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng);
+
+double pIntegrandPS(long lTime, const smc::particle<mChain<double> >& pPos, void* pVoid);
+double pWidthPS(long lTime, void* pVoid);
+double pIntegrandFS(const mChain<double>& dPos, void* pVoid);
+
+///The number of grid elements to either side of the current state for the single state move
+#define GRIDSIZE 12
+///The value of alpha at the specified time
+#define ALPHA(T) (double(T)*double(dSchedule) / double(lIterates))
+///The terminal version of alpha
+#define FTIME (ALPHA(lIterates))
+///The exceedance threshold which we are interested in.
+#define THRESHOLD dThreshold
+///The number of steps in the Markov chain
+#define PATHLENGTH lChainLength
Deleted: pkg/RcppSMC/inst/include/simfunctions.h
===================================================================
--- pkg/RcppSMC/inst/include/simfunctions.h 2012-03-17 02:27:01 UTC (rev 3514)
+++ pkg/RcppSMC/inst/include/simfunctions.h 2012-03-17 09:47:08 UTC (rev 3515)
@@ -1,32 +0,0 @@
-#include "smctc.h"
-
-#include "markovchains/markovchain.h"
-
-extern long lIterates;
-extern long lNumber;
-extern long lChainLength;
-extern double dSchedule;
-extern double dThreshold;
-
-double logDensity(long lTime, const mChain<double> & X);
-
-smc::particle<mChain<double> > fInitialiseMC(smc::rng *pRng);
-long fSelect(long lTime, const smc::particle<mChain<double> > & p, smc::rng *pRng);
-void fMove1(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng);
-void fMove2(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng);
-int fMCMC(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng);
-
-double pIntegrandPS(long lTime, const smc::particle<mChain<double> >& pPos, void* pVoid);
-double pWidthPS(long lTime, void* pVoid);
-double pIntegrandFS(const mChain<double>& dPos, void* pVoid);
-
-///The number of grid elements to either side of the current state for the single state move
-#define GRIDSIZE 12
-///The value of alpha at the specified time
-#define ALPHA(T) (double(T)*double(dSchedule) / double(lIterates))
-///The terminal version of alpha
-#define FTIME (ALPHA(lIterates))
-///The exceedance threshold which we are interested in.
-#define THRESHOLD dThreshold
-///The number of steps in the Markov chain
-#define PATHLENGTH lChainLength
Copied: pkg/RcppSMC/man/rareEventsEx/rareEventsEx.Rd (from rev 3514, pkg/RcppSMC/man/rareEventsEx.Rd)
===================================================================
--- pkg/RcppSMC/man/rareEventsEx/rareEventsEx.Rd (rev 0)
+++ pkg/RcppSMC/man/rareEventsEx/rareEventsEx.Rd 2012-03-17 09:47:08 UTC (rev 3515)
@@ -0,0 +1,46 @@
+\name{rareEventsEx}
+\alias{rareEventsEx}
+\title{Particle Filter Example}
+\description{
+ The \code{rareEventsEx} function provides another example for
+ \pkg{RcppSMC}. It is based on the second example in \code{SMCTC} and
+ the discussion in Section 5.2 of Johansen (2009). It implements an
+ algorithm in Johansen et al (2006) which estimates rare event
+ probabilities for a Gaussian distribution.
+}
+\usage{
+ rareEventsEx(number=100, iterations=10, threshold=5.0, schedule=30.0)
+}
+\arguments{
+ \item{number}{An integer specifying the number of particles.}
+ \item{iterations}{An integer specifying the number of iterations.}
+ \item{threshold}{A numericd describing the rare event threshold.}
+ \item{schedule}{A numericd describing the annealing schedule constant.}
+}
+\value{
+ The function returns a named \code{vector.frame} containing three values.
+}
+\details{
+ The \code{rareEventsEx} function provides the second simple example for
+ \pkg{RcppSMC}. It is based on the \code{rare-events} example in the
+ \code{SMCTC} library, and discussed in the Section 5.2 of his
+ corresponding paper (Johansen, 2009), which in turn refers to Johansen
+ et al (2006) for a deeper discussion.
+}
+\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}
+
+ A. M. Johansen, P. Del Moral, A. Doucet. Sequential Monte Carlo
+ Samplers for Rare Events. Proceedings of the 6th International
+ Workshop on Rare Event Simulation, pages 256--267, Bamberg (Germany).
+}
+\seealso{The SMCTC paper and code at \url{http://www.jstatsoft.org/v30/i06/paper}.}
+\examples{
+\dontrun{
+ res <- rareEventsEx()
+}
+}
+\author{Adam M. Johansen for SMCTC, Dirk Eddelbuettel for the RcppSMC wrapper.}
+\keyword{programming}
Deleted: pkg/RcppSMC/man/rareEventsEx.Rd
===================================================================
--- pkg/RcppSMC/man/rareEventsEx.Rd 2012-03-17 02:27:01 UTC (rev 3514)
+++ pkg/RcppSMC/man/rareEventsEx.Rd 2012-03-17 09:47:08 UTC (rev 3515)
@@ -1,46 +0,0 @@
-\name{rareEventsEx}
-\alias{rareEventsEx}
-\title{Particle Filter Example}
-\description{
- The \code{rareEventsEx} function provides another example for
- \pkg{RcppSMC}. It is based on the second example in \code{SMCTC} and
- the discussion in Section 5.2 of Johansen (2009). It implements an
- algorithm in Johansen et al (2006) which estimates rare event
- probabilities for a Gaussian distribution.
-}
-\usage{
- rareEventsEx(number=100, iterations=10, threshold=5.0, schedule=30.0)
-}
-\arguments{
- \item{number}{An integer specifying the number of particles.}
- \item{iterations}{An integer specifying the number of iterations.}
- \item{threshold}{A numericd describing the rare event threshold.}
- \item{schedule}{A numericd describing the annealing schedule constant.}
-}
-\value{
- The function returns a named \code{vector.frame} containing three values.
-}
-\details{
- The \code{rareEventsEx} function provides the second simple example for
- \pkg{RcppSMC}. It is based on the \code{rare-events} example in the
- \code{SMCTC} library, and discussed in the Section 5.2 of his
- corresponding paper (Johansen, 2009), which in turn refers to Johansen
- et al (2006) for a deeper discussion.
-}
-\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}
-
- A. M. Johansen, P. Del Moral, A. Doucet. Sequential Monte Carlo
- Samplers for Rare Events. Proceedings of the 6th International
- Workshop on Rare Event Simulation, pages 256--267, Bamberg (Germany).
-}
-\seealso{The SMCTC paper and code at \url{http://www.jstatsoft.org/v30/i06/paper}.}
-\examples{
-\dontrun{
- res <- rareEventsEx()
-}
-}
-\author{Adam M. Johansen for SMCTC, Dirk Eddelbuettel for the RcppSMC wrapper.}
-\keyword{programming}
Copied: pkg/RcppSMC/src/rareEvents/rareEvents.cpp (from rev 3514, pkg/RcppSMC/src/rareEvents.cpp)
===================================================================
--- pkg/RcppSMC/src/rareEvents/rareEvents.cpp (rev 0)
+++ pkg/RcppSMC/src/rareEvents/rareEvents.cpp 2012-03-17 09:47:08 UTC (rev 3515)
@@ -0,0 +1,83 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
+//
+// rareEvents.cpp: Rcpp wrapper for SMC library -- second example of Johansen (2009)
+//
+// Copyright (C) 2012 Dirk Eddelbuettel
+//
+// 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 wrap, SMCTC which is
+//
+// Copyright Adam Johansen, 2008.
+//
+// and released under GPL-3, see the copyright headers in inst/include/
+
+#include <Rcpp.h>
+
+#include <iostream>
+#include <cmath>
+#include "simfunctions.h"
+
+///Length of Markov Chain
+long lChainLength = 15;
+///Number of distributions used
+long lIterates;
+///Rare event threshold
+double dThreshold = 5.0;
+///Annealing schedule constant
+double dSchedule = 30.0;
+
+// rareEvents() function callable from R via Rcpp -- which is essentially
+// the same as main() from SMCTC's examples/rare-events/main.cc
+extern "C" SEXP rareEvents(SEXP numberS, SEXP iteratesS, SEXP thresholdS, SEXP scheduleS) {
+
+ long lNumber = Rcpp::as<long>(numberS); // Number of particles
+ lIterates = Rcpp::as<long>(iteratesS); // Number of iterations
+ dThreshold = Rcpp::as<double>(thresholdS); // Rare event threshold
+ dSchedule = Rcpp::as<double>(scheduleS); // Annealing schedule constant
+
+ try{
+ ///An array of move function pointers
+ void (*pfMoves[])(long, smc::particle<mChain<double> > &,smc::rng*) = {fMove1, fMove2};
+ smc::moveset<mChain<double> > Moveset(fInitialiseMC, fSelect, sizeof(pfMoves) / sizeof(pfMoves[0]), pfMoves, fMCMC);
+ smc::sampler<mChain<double> > Sampler(lNumber, SMC_HISTORY_RAM);
+
+ Sampler.SetResampleParams(SMC_RESAMPLE_STRATIFIED,0.5);
+ Sampler.SetMoveSet(Moveset);
+
+ Sampler.Initialise();
+ Sampler.IterateUntil(lIterates);
+
+ ///Estimate the normalising constant of the terminal distribution
+ double zEstimate = Sampler.IntegratePathSampling(pIntegrandPS, pWidthPS, NULL) - log(2.0);
+ ///Estimate the weighting factor for the terminal distribution
+ double wEstimate = Sampler.Integrate(pIntegrandFS, NULL);
+
+ // cout << zEstimate << " " << log(wEstimate) << " " << zEstimate + log(wEstimate) << endl;
+ Rcpp::NumericVector res = Rcpp::NumericVector::create(Rcpp::Named("zEstimate") = zEstimate,
+ Rcpp::Named("log(wEstimate)") = log(wEstimate),
+ Rcpp::Named("sum") = zEstimate + log(wEstimate));
+ return res;
+ }
+ 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
+
+}
+
Copied: pkg/RcppSMC/src/rareEvents/simfunctions.cpp (from rev 3514, pkg/RcppSMC/src/simfunctions.cpp)
===================================================================
--- pkg/RcppSMC/src/rareEvents/simfunctions.cpp (rev 0)
+++ pkg/RcppSMC/src/rareEvents/simfunctions.cpp 2012-03-17 09:47:08 UTC (rev 3515)
@@ -0,0 +1,171 @@
+#include <iostream>
+#include <cmath>
+//#include <gsl/gsl_randist.h>
+#include <Rcpp.h>
+
+#include "smctc.h"
+#include "simfunctions.h"
+
+using namespace std;
+
+///The function corresponding to the log posterior density at specified time and position
+
+/// \param lTime The current time (i.e. the index of the current distribution)
+/// \param X The state to consider **/
+double logDensity(long lTime, const mChain<double> & X)
+{
+ double lp;
+
+ mElement<double> *x = X.GetElement(0);
+ mElement<double> *y = x->pNext;
+ //Begin with the density exluding the effect of the potential
+ //lp = log(gsl_ran_ugaussian_pdf(x->value));
+ Rcpp::NumericVector z = Rcpp::dnorm( Rcpp::NumericVector(x->value), 0.0, 1.0);
+ lp = log( Rcpp::as<double>( z ));
+
+ while(y) {
+ //lp += log(gsl_ran_ugaussian_pdf(y->value - x->value));
+ z = Rcpp::dnorm( Rcpp::NumericVector(y->value - x->value), 0.0, 1.0);
+ lp += log( Rcpp::as<double>(z));
+ x = y;
+ y = x->pNext;
+ }
+
+ //Now include the effect of the multiplicative potential function
+ lp -= log(1.0 + exp(-(ALPHA(lTime) * (x->value - THRESHOLD) )));
+ return lp;
+}
+
+///A function to initialise double type markov chain-valued particles
+/// \param pRng A pointer to the random number generator which is to be used
+smc::particle<mChain<double> > fInitialiseMC(smc::rng *pRng)
+{
+ // Create a Markov chain with the appropriate initialisation and then assign that to the particles.
+ mChain<double> Mc;
+
+ double x = 0;
+ for(int i = 0; i < PATHLENGTH; i++) {
+ x += pRng->NormalS();
+ Mc.AppendElement(x);
+ }
+
+ return smc::particle<mChain<double> >(Mc,0);
+}
+
+///A function to select a move randomly
+/// \param lTime The current evolution time of the system
+/// \param p The current position of the particle which is to be moved
+/// \param pRng A pointer to the random number generator which is to be used
+long fSelect(long lTime, const smc::particle<mChain<double> > & p, smc::rng *pRng)
+{
+ return 0;
+}
+
+void fMove1(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng)
+{
+ // The distance between points in the random grid.
+ static double delta = 0.025;
+ static double gridweight[2*GRIDSIZE+1], gridws = 0;
+ static mChain<double> NewPos[2*GRIDSIZE+1];
+ static mChain<double> OldPos[2*GRIDSIZE+1];
+
+ // First select a new position from a grid centred on the old position, weighting the possible choises by the
+ // posterior probability of the resulting states.
+ gridws = 0;
+ for(int i = 0; i < 2*GRIDSIZE+1; i++) {
+ NewPos[i] = pFrom.GetValue() + ((double)(i - GRIDSIZE))*delta;
+ gridweight[i] = exp(logDensity(lTime,NewPos[i]));
+ gridws = gridws + gridweight[i];
+ }
+
+ double dRUnif = pRng->Uniform(0,gridws);
+ long j = -1;
+
+ while(dRUnif > 0 && j <= 2*GRIDSIZE) {
+ j++;
+ dRUnif -= gridweight[j];
+ }
+
+ pFrom.SetValue(NewPos[j]);
+
+ // Now calculate the weight change which the particle suffers as a result
+ double logInc = log(gridweight[j]), Inc = 0;
+
+ for(int i = 0; i < 2*GRIDSIZE+1; i++) {
+ OldPos[i] = pFrom.GetValue() - ((double)(i - GRIDSIZE))*delta;
+ gridws = 0;
+ for(int k = 0; k < 2*GRIDSIZE+1; k++) {
+ NewPos[k] = OldPos[i] + ((double)(k-GRIDSIZE))*delta;
+ gridweight[k] = exp(logDensity(lTime, NewPos[k]));
+ gridws += gridweight[k];
+ }
+ Inc += exp(logDensity(lTime-1, OldPos[i])) * exp(logDensity(lTime, pFrom.GetValue())) / gridws;
+ }
+ logInc -= log(Inc);
+
+ pFrom.SetLogWeight(pFrom.GetLogWeight() + logInc);
+
+ for(int i = 0; i < 2*GRIDSIZE+1; i++)
+ {
+ NewPos[i].Empty();
+ OldPos[i].Empty();
+ }
+
+ return;
+}
+///Another move function
+void fMove2(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng)
+{
+ pFrom.SetLogWeight(pFrom.GetLogWeight() + logDensity(lTime,pFrom.GetValue()) - logDensity(lTime-1,pFrom.GetValue()));
+}
+
+///An MCMC step suitable for introducing sample diversity
+int fMCMC(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng)
+{
+ static smc::particle<mChain<double> > pTo;
+
+ mChain<double> * pMC = new mChain<double>;
+
+ for(int i = 0; i < pFrom.GetValue().GetLength(); i++)
+ pMC->AppendElement(pFrom.GetValue().GetElement(i)->value + pRng->Normal(0, 0.5));
+ pTo.SetValue(*pMC);
+ pTo.SetLogWeight(pFrom.GetLogWeight());
+
+ delete pMC;
+
+ double alpha = exp(logDensity(lTime,pTo.GetValue()) - logDensity(lTime,pFrom.GetValue()));
+ if(alpha < 1)
+ if (pRng->UniformS() > alpha) {
+ return false;
+ }
+
+ pFrom = pTo;
+ return true;
+}
+
+///A function to be integrated in the path sampling step.
+double pIntegrandPS(long lTime, const smc::particle<mChain<double> >& pPos, void* pVoid)
+{
+ double dPos = pPos.GetValue().GetTerminal()->value;
+ return (dPos - THRESHOLD) / (1.0 + exp(ALPHA(lTime) * (dPos - THRESHOLD)));
+}
+
+///A function which gives the width distribution for the path sampling step.
+double pWidthPS(long lTime, void* pVoid)
+{
+ if(lTime > 1 && lTime < lIterates)
+ return ((0.5)*double(ALPHA(lTime+1.0)-ALPHA(lTime-1.0)));
+ else
+ return((0.5)*double(ALPHA(lTime+1.0)-ALPHA(lTime)) +(ALPHA(1)-0.0));
+}
+
+//The final state weighting function -- how likely is a random path from this distribution to hit the rare set...
+double pIntegrandFS(const mChain<double>& dPos, void* pVoid)
+{
+ if(dPos.GetTerminal()->value > THRESHOLD) {
+ return (1.0 + exp(-FTIME*(dPos.GetTerminal()->value-THRESHOLD)));
+ }
+ else
+ return 0;
+}
+
Deleted: pkg/RcppSMC/src/rareEvents.cpp
===================================================================
--- pkg/RcppSMC/src/rareEvents.cpp 2012-03-17 02:27:01 UTC (rev 3514)
+++ pkg/RcppSMC/src/rareEvents.cpp 2012-03-17 09:47:08 UTC (rev 3515)
@@ -1,83 +0,0 @@
-// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
-//
-// rareEvents.cpp: Rcpp wrapper for SMC library -- second example of Johansen (2009)
-//
-// Copyright (C) 2012 Dirk Eddelbuettel
-//
-// 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 wrap, SMCTC which is
-//
-// Copyright Adam Johansen, 2008.
-//
-// and released under GPL-3, see the copyright headers in inst/include/
-
-#include <Rcpp.h>
-
-#include <iostream>
-#include <cmath>
-#include "simfunctions.h"
-
-///Length of Markov Chain
-long lChainLength = 15;
-///Number of distributions used
-long lIterates;
-///Rare event threshold
-double dThreshold = 5.0;
-///Annealing schedule constant
-double dSchedule = 30.0;
-
-// rareEvents() function callable from R via Rcpp -- which is essentially
-// the same as main() from SMCTC's examples/rare-events/main.cc
-extern "C" SEXP rareEvents(SEXP numberS, SEXP iteratesS, SEXP thresholdS, SEXP scheduleS) {
-
- long lNumber = Rcpp::as<long>(numberS); // Number of particles
- lIterates = Rcpp::as<long>(iteratesS); // Number of iterations
- dThreshold = Rcpp::as<double>(thresholdS); // Rare event threshold
- dSchedule = Rcpp::as<double>(scheduleS); // Annealing schedule constant
-
- try{
- ///An array of move function pointers
- void (*pfMoves[])(long, smc::particle<mChain<double> > &,smc::rng*) = {fMove1, fMove2};
- smc::moveset<mChain<double> > Moveset(fInitialiseMC, fSelect, sizeof(pfMoves) / sizeof(pfMoves[0]), pfMoves, fMCMC);
- smc::sampler<mChain<double> > Sampler(lNumber, SMC_HISTORY_RAM);
-
- Sampler.SetResampleParams(SMC_RESAMPLE_STRATIFIED,0.5);
- Sampler.SetMoveSet(Moveset);
-
- Sampler.Initialise();
- Sampler.IterateUntil(lIterates);
-
- ///Estimate the normalising constant of the terminal distribution
- double zEstimate = Sampler.IntegratePathSampling(pIntegrandPS, pWidthPS, NULL) - log(2.0);
- ///Estimate the weighting factor for the terminal distribution
- double wEstimate = Sampler.Integrate(pIntegrandFS, NULL);
-
- // cout << zEstimate << " " << log(wEstimate) << " " << zEstimate + log(wEstimate) << endl;
- Rcpp::NumericVector res = Rcpp::NumericVector::create(Rcpp::Named("zEstimate") = zEstimate,
- Rcpp::Named("log(wEstimate)") = log(wEstimate),
- Rcpp::Named("sum") = zEstimate + log(wEstimate));
- return res;
- }
- 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
-
-}
-
Deleted: pkg/RcppSMC/src/simfunctions.cpp
===================================================================
--- pkg/RcppSMC/src/simfunctions.cpp 2012-03-17 02:27:01 UTC (rev 3514)
+++ pkg/RcppSMC/src/simfunctions.cpp 2012-03-17 09:47:08 UTC (rev 3515)
@@ -1,171 +0,0 @@
-#include <iostream>
-#include <cmath>
-//#include <gsl/gsl_randist.h>
-#include <Rcpp.h>
-
-#include "smctc.h"
-#include "simfunctions.h"
-
-using namespace std;
-
-///The function corresponding to the log posterior density at specified time and position
-
-/// \param lTime The current time (i.e. the index of the current distribution)
-/// \param X The state to consider **/
-double logDensity(long lTime, const mChain<double> & X)
-{
- double lp;
-
- mElement<double> *x = X.GetElement(0);
- mElement<double> *y = x->pNext;
- //Begin with the density exluding the effect of the potential
- //lp = log(gsl_ran_ugaussian_pdf(x->value));
- Rcpp::NumericVector z = Rcpp::dnorm( Rcpp::NumericVector(x->value), 0.0, 1.0);
- lp = log( Rcpp::as<double>( z ));
-
- while(y) {
- //lp += log(gsl_ran_ugaussian_pdf(y->value - x->value));
- z = Rcpp::dnorm( Rcpp::NumericVector(y->value - x->value), 0.0, 1.0);
- lp += log( Rcpp::as<double>(z));
- x = y;
- y = x->pNext;
- }
-
- //Now include the effect of the multiplicative potential function
- lp -= log(1.0 + exp(-(ALPHA(lTime) * (x->value - THRESHOLD) )));
- return lp;
-}
-
-///A function to initialise double type markov chain-valued particles
-/// \param pRng A pointer to the random number generator which is to be used
-smc::particle<mChain<double> > fInitialiseMC(smc::rng *pRng)
-{
- // Create a Markov chain with the appropriate initialisation and then assign that to the particles.
- mChain<double> Mc;
-
- double x = 0;
- for(int i = 0; i < PATHLENGTH; i++) {
- x += pRng->NormalS();
- Mc.AppendElement(x);
- }
-
- return smc::particle<mChain<double> >(Mc,0);
-}
-
-///A function to select a move randomly
-/// \param lTime The current evolution time of the system
-/// \param p The current position of the particle which is to be moved
-/// \param pRng A pointer to the random number generator which is to be used
-long fSelect(long lTime, const smc::particle<mChain<double> > & p, smc::rng *pRng)
-{
- return 0;
-}
-
-void fMove1(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng)
-{
- // The distance between points in the random grid.
- static double delta = 0.025;
- static double gridweight[2*GRIDSIZE+1], gridws = 0;
- static mChain<double> NewPos[2*GRIDSIZE+1];
- static mChain<double> OldPos[2*GRIDSIZE+1];
-
- // First select a new position from a grid centred on the old position, weighting the possible choises by the
- // posterior probability of the resulting states.
- gridws = 0;
- for(int i = 0; i < 2*GRIDSIZE+1; i++) {
- NewPos[i] = pFrom.GetValue() + ((double)(i - GRIDSIZE))*delta;
- gridweight[i] = exp(logDensity(lTime,NewPos[i]));
- gridws = gridws + gridweight[i];
- }
-
- double dRUnif = pRng->Uniform(0,gridws);
- long j = -1;
-
- while(dRUnif > 0 && j <= 2*GRIDSIZE) {
- j++;
- dRUnif -= gridweight[j];
- }
-
- pFrom.SetValue(NewPos[j]);
-
- // Now calculate the weight change which the particle suffers as a result
- double logInc = log(gridweight[j]), Inc = 0;
-
- for(int i = 0; i < 2*GRIDSIZE+1; i++) {
- OldPos[i] = pFrom.GetValue() - ((double)(i - GRIDSIZE))*delta;
- gridws = 0;
- for(int k = 0; k < 2*GRIDSIZE+1; k++) {
- NewPos[k] = OldPos[i] + ((double)(k-GRIDSIZE))*delta;
- gridweight[k] = exp(logDensity(lTime, NewPos[k]));
- gridws += gridweight[k];
- }
- Inc += exp(logDensity(lTime-1, OldPos[i])) * exp(logDensity(lTime, pFrom.GetValue())) / gridws;
- }
- logInc -= log(Inc);
-
- pFrom.SetLogWeight(pFrom.GetLogWeight() + logInc);
-
- for(int i = 0; i < 2*GRIDSIZE+1; i++)
- {
- NewPos[i].Empty();
- OldPos[i].Empty();
- }
-
- return;
-}
-///Another move function
-void fMove2(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng)
-{
- pFrom.SetLogWeight(pFrom.GetLogWeight() + logDensity(lTime,pFrom.GetValue()) - logDensity(lTime-1,pFrom.GetValue()));
-}
-
-///An MCMC step suitable for introducing sample diversity
-int fMCMC(long lTime, smc::particle<mChain<double> > & pFrom, smc::rng *pRng)
-{
- static smc::particle<mChain<double> > pTo;
-
- mChain<double> * pMC = new mChain<double>;
-
- for(int i = 0; i < pFrom.GetValue().GetLength(); i++)
- pMC->AppendElement(pFrom.GetValue().GetElement(i)->value + pRng->Normal(0, 0.5));
- pTo.SetValue(*pMC);
- pTo.SetLogWeight(pFrom.GetLogWeight());
-
- delete pMC;
-
- double alpha = exp(logDensity(lTime,pTo.GetValue()) - logDensity(lTime,pFrom.GetValue()));
- if(alpha < 1)
- if (pRng->UniformS() > alpha) {
- return false;
- }
-
- pFrom = pTo;
- return true;
-}
-
-///A function to be integrated in the path sampling step.
-double pIntegrandPS(long lTime, const smc::particle<mChain<double> >& pPos, void* pVoid)
-{
- double dPos = pPos.GetValue().GetTerminal()->value;
- return (dPos - THRESHOLD) / (1.0 + exp(ALPHA(lTime) * (dPos - THRESHOLD)));
-}
-
-///A function which gives the width distribution for the path sampling step.
-double pWidthPS(long lTime, void* pVoid)
-{
- if(lTime > 1 && lTime < lIterates)
- return ((0.5)*double(ALPHA(lTime+1.0)-ALPHA(lTime-1.0)));
- else
- return((0.5)*double(ALPHA(lTime+1.0)-ALPHA(lTime)) +(ALPHA(1)-0.0));
-}
-
-//The final state weighting function -- how likely is a random path from this distribution to hit the rare set...
-double pIntegrandFS(const mChain<double>& dPos, void* pVoid)
-{
- if(dPos.GetTerminal()->value > THRESHOLD) {
- return (1.0 + exp(-FTIME*(dPos.GetTerminal()->value-THRESHOLD)));
- }
- else
- return 0;
-}
-
More information about the Rcpp-commits
mailing list