[Rcpp-commits] r3441 - in pkg/RcppSMC: . R man src src/markovchains
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jan 16 00:16:13 CET 2012
Author: edd
Date: 2012-01-16 00:16:13 +0100 (Mon, 16 Jan 2012)
New Revision: 3441
Added:
pkg/RcppSMC/R/rareEventsEx.R
pkg/RcppSMC/man/rareEventsEx.Rd
pkg/RcppSMC/src/markovchains/
pkg/RcppSMC/src/markovchains/markovchain.h
pkg/RcppSMC/src/rareEvents.cpp
pkg/RcppSMC/src/simfunctions.cc
pkg/RcppSMC/src/simfunctions.hh
Modified:
pkg/RcppSMC/ChangeLog
pkg/RcppSMC/NAMESPACE
pkg/RcppSMC/src/pf.cpp
Log:
added second example from Johansen (2009)
Modified: pkg/RcppSMC/ChangeLog
===================================================================
--- pkg/RcppSMC/ChangeLog 2012-01-15 10:42:22 UTC (rev 3440)
+++ pkg/RcppSMC/ChangeLog 2012-01-15 23:16:13 UTC (rev 3441)
@@ -1,3 +1,14 @@
+2012-01-15 Dirk Eddelbuettel <edd at debian.org>
+
+ * src/rareEventsEx.cpp: Adapted main() from rare-events/main.cc in
+ the SMCTC sources
+ * src/simfunctions.hh: Also imported from rare-events/, and renamed
+ fInitialize to fInitializMC as the library already fInitialize
+ * src/simfunctions.hh: Dito
+ * src/markovchains/markovchain.h: Also imported from rare-events/
+ * R/rareEvents.R: Added minimal wrapper function
+ * man/rareEvents.Rd: Added minimal manual page
+
2012-01-11 Dirk Eddelbuettel <edd at debian.org>
* Initial commit to R-Forge; package is still very bare but at least
Modified: pkg/RcppSMC/NAMESPACE
===================================================================
--- pkg/RcppSMC/NAMESPACE 2012-01-15 10:42:22 UTC (rev 3440)
+++ pkg/RcppSMC/NAMESPACE 2012-01-15 23:16:13 UTC (rev 3441)
@@ -1,4 +1,5 @@
useDynLib(RcppSMC)
-export(pfEx)
+export(pfEx,
+ rareEventsEx)
Added: pkg/RcppSMC/R/rareEventsEx.R
===================================================================
--- pkg/RcppSMC/R/rareEventsEx.R (rev 0)
+++ pkg/RcppSMC/R/rareEventsEx.R 2012-01-15 23:16:13 UTC (rev 3441)
@@ -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)
+}
+
Added: pkg/RcppSMC/man/rareEventsEx.Rd
===================================================================
--- pkg/RcppSMC/man/rareEventsEx.Rd (rev 0)
+++ pkg/RcppSMC/man/rareEventsEx.Rd 2012-01-15 23:16:13 UTC (rev 3441)
@@ -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}
Added: pkg/RcppSMC/src/markovchains/markovchain.h
===================================================================
--- pkg/RcppSMC/src/markovchains/markovchain.h (rev 0)
+++ pkg/RcppSMC/src/markovchains/markovchain.h 2012-01-15 23:16:13 UTC (rev 3441)
@@ -0,0 +1,593 @@
+#ifndef __MARKOVCHAIN_H
+#define __MARKOVCHAIN_H
+
+#include <cstdlib>
+#include <iostream>
+
+using namespace std;
+
+//!A template class for the elements of a Markov chain.
+template <class Space> class mElement
+{
+ public:
+ //! Pointer to the previous element in the chain.
+ mElement *pPrev;
+ //! Pointer to the next element in the chain.
+ mElement *pNext;
+ //! Value of this element of the chain.
+ Space value;
+
+ mElement();
+ mElement(Space);
+ mElement(mElement<Space>*,Space,mElement<Space>*);
+
+ mElement<Space> operator=(mElement<Space> & mE) { this->value = mE.value; }
+ mElement<Space> operator=(mElement<Space> mE) { this->value = mE.value; }
+};
+
+
+template <class Space>
+mElement<Space>::mElement()
+{
+ pPrev = NULL;
+ pNext = NULL;
+}
+
+template <class Space>
+mElement<Space>::mElement(Space sInit)
+{
+ value = sInit;
+
+ pPrev = NULL;
+ pNext = NULL;
+}
+
+template <class Space>
+mElement<Space>::mElement(mElement<Space>* pPInit, Space sInit, mElement<Space>* pNInit)
+{
+ value = sInit;
+
+ pPrev = pPInit;
+ pNext = pNInit;
+}
+
+//!A template class for a Markov Chain itself.
+template <class Space> class mChain
+{
+ private:
+ long nLength;
+ mElement<Space> *pRoot;
+ mElement<Space> *pLeaf;
+
+ Space (*pfKernel) (Space);
+ public:
+ mChain();
+ mChain(long nInit, Space* pElements);
+ mChain(mChain<Space> *pMci);
+ mChain(const mChain<Space> &mci);
+
+ ~mChain();
+
+ mChain<Space> operator= (const mChain<Space>&);
+ mChain<Space> operator+ (Space const &) const;
+ mChain<Space>& operator+=(Space const &);
+ mChain<Space> operator- (Space const &) const;
+ mChain<Space>& operator-=(Space const &);
+ const Space & AddToElement(long nIndex, const Space &);
+ int AppendElement(const Space & );
+ int CheckSanity(void) const;
+ int DeleteElement(long);
+ int DeleteRange(long,long);
+ void Empty(void);
+ int KernelElement(void);
+ int KernelRange(long);
+ int InsertElement(long, Space);
+ int PrependElement(Space);
+
+ mElement<Space> *GetElement(long ind) const {if(ind > nLength) return NULL; mElement<Space> *pCount = pRoot; for(int i=0; i < ind; i++) pCount = pCount->pNext; return pCount; }
+ mElement<Space> const *GetTerminal(void) const {return pLeaf;}
+ long GetLength(void) const { return nLength; }
+
+ Space (*GetKernel(void))(Space) { return pfKernel; }
+ void SetKernel(Space (*pfKernelNew) (Space)) { pfKernel = pfKernelNew; }
+};
+
+//! Constructor for an empty Markov Chain.
+template <class Space>
+mChain<Space>::mChain()
+{
+ nLength = 0;
+ pRoot = NULL;
+ pLeaf = NULL;
+
+ pfKernel = NULL;
+
+}
+
+//! Constructor for an initialised Markov Chain.
+template <class Space>
+mChain<Space>::mChain(long nInit, Space* pElements)
+{
+ nLength = nInit;
+
+ if(nInit) {
+ pRoot = new mElement<Space>(pElements[0]);
+
+ mElement<Space> *pCurrent = pRoot;
+
+ for(int i = 1; i < nInit; i++) {
+ pCurrent->pNext = new mElement<Space>(pCurrent, pElements[i], NULL);
+ pCurrent = pCurrent->pNext;
+ }
+
+ pLeaf = pCurrent;
+ }
+ else {
+ pRoot = NULL;
+ pLeaf = NULL;
+ }
+
+ pfKernel = NULL;
+}
+
+//! Constructor for an initialisd Markov Chain
+template <class Space>
+mChain<Space>::mChain(mChain<Space> *pMci)
+{
+ nLength = pMci->mLength;
+ if(pMci->nLength == 0) {
+ // nLength = 0;
+ pRoot = pLeaf = NULL;
+ }
+ else if(pMci->nLength == 1) {
+ pRoot = new mElement<Space>(NULL, pMci->GetElement(0)->value,NULL);
+ pLeaf = pRoot;
+ // nLength = 1;
+ }
+ else {
+ mElement<Space> * mEN = pMci->GetElement(0);
+ pRoot = new mElement<Space>(NULL,mEN->value,NULL);
+ mElement<Space> * mCurrent = pRoot;
+ for(int i = 1; i < nLength; i++) {
+ mCurrent->pNext = new mElement<Space>(mCurrent, mEN->value, NULL);
+ mEN = mEN->pNext;
+ mCurrent = mCurrent->pNext;
+ }
+ pLeaf = mCurrent;
+ }
+
+}
+
+template <class Space>
+mChain<Space>::mChain(const mChain<Space> &Mci)
+{
+ /*
+ nLength = 0;
+ pRoot = pLeaf = NULL;
+
+ for(int i = 0; i < Mci.nLength; i++)
+ AppendElement(Mci.GetElement(i)->value);
+ */
+
+ nLength = Mci.nLength;
+ if(Mci.nLength == 0) {
+ pRoot = pLeaf = NULL;
+ }
+ else {
+ if(Mci.nLength == 1) {
+ pRoot = new mElement<Space>(NULL, Mci.GetElement(0)->value,NULL);
+ pLeaf = pRoot;
+ }
+ else {
+ mElement<Space> * mEN = Mci.GetElement(0);
+ pRoot = new mElement<Space>(NULL,mEN->value,NULL);
+ mElement<Space> * mCurrent = pRoot;
+ for(int i = 1; i < nLength; i++) {
+ mCurrent->pNext = new mElement<Space>(mCurrent, mEN->pNext->value, NULL);
+ mEN = mEN->pNext;
+ mCurrent = mCurrent->pNext;
+ }
+ pLeaf = mCurrent;
+ }
+ }
+
+ pfKernel = Mci.pfKernel;
+}
+
+template <class Space>
+mChain<Space>::~mChain()
+{
+ //It's wise to have something here or things leak somewhat....
+
+ if(nLength > 0) {
+ mElement<Space>* pCurrent = pLeaf;
+ while(pCurrent->pPrev) {
+ pCurrent = pCurrent->pPrev;
+ pCurrent->pNext->pPrev = NULL;
+ delete pCurrent->pNext;
+ pCurrent->pNext=NULL;
+ }
+ delete pCurrent;
+ }
+ pRoot=pLeaf=NULL;
+ nLength = 0;
+}
+
+/// Add the specified value to the specified element of the Markov Chain
+
+template <class Space>
+const Space & mChain<Space>::AddToElement(long nElement, const Space & sAdd)
+{
+ mElement<Space>* pCurrent = pRoot;
+
+ for(int n = 0; n < nElement; n++)
+ {
+ pCurrent=pCurrent->pNext;
+ if(!pCurrent)
+ return sAdd;
+ }
+
+ return pCurrent->value = pCurrent->value + sAdd;
+}
+
+//! Add a specified element to the end of a Markov Chain.
+template <class Space>
+int inline mChain<Space>::AppendElement(const Space &snew)
+{
+ if(nLength) {
+ pLeaf->pNext = new mElement<Space>(pLeaf,snew,NULL);
+ pLeaf = pLeaf->pNext;
+
+ nLength++;
+
+ return 0;
+ }
+ else {
+ pRoot = new mElement<Space>(NULL,snew,NULL);
+ pLeaf = pRoot;
+ nLength = 1;
+ return 0;
+ }
+}
+
+//! Check that the doubly-linked list is sane
+template <class Space>
+int mChain<Space>::CheckSanity(void) const
+{
+ int nLengthCheck = 0, nLengthCheck2 = 0;
+ int nLinkFails = 0;
+
+ mElement<Space> *pHead = pRoot;
+
+ while(pHead) {
+ nLengthCheck++;
+
+ if(pHead->pNext) {
+ if(pHead->pNext->pPrev != pHead)
+ nLinkFails++;
+ }
+
+ pHead=pHead->pNext;
+ }
+
+ pHead = pLeaf;
+ while(pHead) {
+ nLengthCheck2++;
+
+ if(pHead->pPrev) {
+ if(pHead->pPrev->pNext != pHead)
+ nLinkFails++;
+ }
+
+ pHead=pHead->pPrev;
+ }
+
+ if(nLengthCheck != nLength || nLengthCheck2 != nLength)
+ cerr << "LengthCheck " << nLengthCheck << "," << nLengthCheck2 << "(" << nLength << ")" << endl;
+ if(nLinkFails)
+ cerr << "Link Failure count..." << nLinkFails << endl;
+
+ if(nLengthCheck != nLength)
+ return -1;
+ return nLinkFails;
+}
+
+//! Remove an element from a specified position in a Markov Chain.
+template <class Space>
+int mChain<Space>::DeleteElement(long nPos)
+{
+ if(nPos < 0 || nPos > nLength)
+ return 1;
+
+ if(nPos == 0) {
+ mElement<Space> *pRootN = pRoot->pNext;
+ delete pRoot;
+ pRoot=pRootN;
+ if(pRoot)
+ pRoot->pPrev = NULL;
+
+ nLength--;
+
+ if(nLength < 2)
+ pLeaf = pRoot;
+
+ return 0;
+ }
+ if(nPos == nLength) {
+ mElement<Space> *pLeafN = pLeaf->pPrev;
+ delete pLeaf;
+ pLeaf=pLeafN;
+ if(pLeaf)
+ pLeaf->pNext = NULL;
+
+ nLength--;
+
+ if(nLength < 2)
+ pRoot = pLeaf;
+
+ return 0;
+ }
+
+ mElement<Space> *pPoint = pRoot, *pPointN, *pPointP;
+ for(int i = 1; i < nPos; i++)
+ pPoint = pPoint->pNext;
+ pPointP = pPoint->pPrev;
+ pPointN = pPoint->pNext;
+
+ delete pPoint;
+ pPointP->pNext = pPointN;
+ pPointN->pPrev = pPointP;
+
+ nLength--;
+
+ return 0;
+}
+
+//! Delete a contiguous region from within a Markov Chain.
+template <class Space>
+int mChain<Space>::DeleteRange(long nPos1, long nPos2)
+{
+ for(int i = nPos1; i < nPos2; i++)
+ if(DeleteElement(nPos1))
+ return 1;
+
+ return 0;
+
+
+ //This should be replaced wite an efficient implementation at some point!
+}
+
+//! Remove all of the elements from a Markov chain, leaving only an empty container
+template <class Space>
+void inline mChain<Space>::Empty(void)
+{
+ if(nLength > 0) {
+ mElement<Space>* pCurrent = pLeaf;
+ while(pCurrent->pPrev) {
+ pCurrent = pCurrent->pPrev;
+ pCurrent->pNext->pPrev = NULL;
+ delete pCurrent->pNext;
+ pCurrent->pNext=NULL;
+ }
+
+ delete pCurrent;
+ }
+ pRoot=pLeaf=NULL;
+ nLength = 0;
+}
+
+//! Insert a specified element at a specified position in a Markov Chain.
+template <class Space>
+int mChain<Space>::InsertElement(long nPos, Space snew)
+{
+ //If the insertion position is zero then insert the new element at the root of the chain
+ if(nPos == 0) {
+ return PrependElement(snew);
+ }
+ if(nPos == nLength)
+ return AppendElement(snew);
+
+ if(nPos > 0 && nPos < nLength) {
+ mElement<Space> *pPoint = pRoot, *pPointN;
+ for(int i = 1; i < nPos; i++)
+ pPoint = pPoint->pNext;
+ pPointN = pPoint->pNext;
+ pPoint->pNext = new mElement<Space>(pPoint,snew,pPointN);
+ pPointN->pPrev= pPoint->pNext;
+
+ nLength++;
+ return 0;
+ }
+
+ //If the insertion position does not lie between 0 and the length of the chain then we have a problem
+ return 1;
+}
+
+//! Extend the chain by a single element obtained from the associated Markov Kernel
+template <class Space>
+int mChain<Space>::KernelElement(void)
+{
+ if(nLength) {
+ pLeaf->pNext = new mElement<Space>(pLeaf, (*pfKernel)(pLeaf->value),NULL);
+ pLeaf = pLeaf->pNext;
+
+ nLength++;
+ return 0;
+ }
+ else {
+ return 1; //We can't very well extend a non-existant chain according to a kernel
+ }
+}
+
+//! Extend the chain by a specified number of elements obtained by iterative application of the associated Markov Kernel.
+template <class Space>
+int mChain<Space>::KernelRange(long n)
+{
+ for(int i = 0; i < n; i++) {
+ if(KernelElement())
+ return 1+i;
+ }
+
+ return 0;
+}
+
+//! Insert a specified element at the beginning of a Markov Chain.
+template <class Space>
+int mChain<Space>::PrependElement(Space snew)
+{
+ if(nLength) {
+ pRoot->pPrev = new mElement<Space>(NULL,snew,pRoot);
+ pRoot=pRoot->pPrev;
+
+ nLength++;
+ return 0;
+ }
+ else {
+ pRoot = new mElement<Space>(NULL,snew,NULL);
+ pLeaf = pRoot;
+ nLength = 1;
+ return 0;
+ }
+}
+
+
+//Overloaded operators for use with these template classes
+
+//! Assigning a Markov Chain performs a deep copy
+/*
+template <class Space>
+mChain<Space> & mChain<Space>::operator= (mChain<Space> & mci)
+{
+ nLength = mci.nLength;
+
+ // mChain<Space> *that = new mChain<Space>(&mci);
+ //
+ mElement<Space> *pi, *po;
+
+ pi = mci.pRoot;
+
+ Empty();
+
+ pRoot = new mElement<Space>(NULL,pi->value,NULL);
+ po = pRoot;
+ // po->value = pi->value;
+
+ while(pi->pNext) {
+ po->pNext = new mElement<Space>(po,pi->pNext->value,NULL);
+ // po->pNext->value = pi->pNext->value;
+ // po->pNext->pPrev = po;
+ po=po->pNext;
+ pi=pi->pNext;
+ }
+ pLeaf = pi;
+
+ return (*this);
+
+ //return *that;
+}
+*/
+ //! Assigning a Markov Chain performs a deep copy
+template <class Space>
+mChain<Space> mChain<Space>::operator= (const mChain<Space> &mci)
+{
+ // First, get rid of the old stuff
+ Empty();
+
+ // Now assign some new stuff
+
+ nLength = mci.nLength;
+
+ if(!nLength)
+ return *this;
+
+ // mChain<Space> *that = new mChain<Space>(&mci);
+ //
+ mElement<Space> *pi, *po;
+
+ pi = mci.pRoot;
+
+ pRoot = new mElement<Space>(NULL,pi->value,NULL);
+ po = pRoot;
+
+ while(pi->pNext) {
+ pi=pi->pNext;
+ po->pNext = new mElement<Space>(po,pi->value,NULL);
+ po=po->pNext;
+ }
+ pLeaf = po;
+
+ return (*this);
+
+ //return *that;
+}
+
+// Provide shortcuts for shifting the entire chain
+template <class Space>
+mChain<Space> mChain<Space>::operator+ (Space const & sInc) const
+{
+ static mChain<Space> that; // = new mChain<Space>;
+
+ that = *this;
+
+ mElement<Space> *pCurrent;
+ pCurrent = that.pRoot;
+ do {
+ pCurrent->value = pCurrent->value + sInc;
+ pCurrent = pCurrent->pNext;
+ } while (pCurrent);
+
+ return that;
+}
+
+template <class Space>
+mChain<Space>& mChain<Space>::operator+=(Space const & sInc)
+{
+ *this = *this + sInc;
+
+ return *this;
+}
+
+template <class Space>
+mChain<Space> mChain<Space>::operator- (Space const & sInc) const
+{
+ static mChain<Space> that; // = new mChain<Space>;
+
+ that = *this;
+
+ mElement<Space> *pCurrent;
+ pCurrent = that.pRoot;
+ do {
+ pCurrent->value = pCurrent->value - sInc;
+ pCurrent = pCurrent->pNext;
+ } while (pCurrent);
+
+ return that;
+}
+
+template <class Space>
+mChain<Space>& mChain<Space>::operator-=(Space const & sInc)
+{
+ *this = *this - sInc;
+
+ return *this;
+}
+//Overloaded stream operators for use with these template classes
+//! Allow the use of output streams with Markov Chains.
+template <class Space>
+ostream & operator<< (ostream& os, mChain<Space>& mc)
+{
+ if(!mc.GetLength())
+ return os;
+
+ mElement<Space> *me = mc.GetElement(0);
+
+ os << "[" << me->value;
+ while(me->pNext) {
+ me = me->pNext;
+ os << "," << me->value;
+ }
+ os << "]";
+ return os;
+}
+
+#endif
Modified: pkg/RcppSMC/src/pf.cpp
===================================================================
--- pkg/RcppSMC/src/pf.cpp 2012-01-15 10:42:22 UTC (rev 3440)
+++ pkg/RcppSMC/src/pf.cpp 2012-01-15 23:16:13 UTC (rev 3441)
@@ -1,6 +1,6 @@
// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
//
-// pf.cpp: Rcpp wrapper for SMC library -- first example
+// pf.cpp: Rcpp wrapper for SMC library -- first example of Johansen (2009)
//
// Copyright (C) 2012 Dirk Eddelbuettel
//
@@ -20,7 +20,7 @@
// 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 and wrap SMCTC which is
+// RcppSMC builds on top of, and wrap, SMCTC which is
//
// Copyright Adam Johansen, 2008.
//
Added: pkg/RcppSMC/src/rareEvents.cpp
===================================================================
--- pkg/RcppSMC/src/rareEvents.cpp (rev 0)
+++ pkg/RcppSMC/src/rareEvents.cpp 2012-01-15 23:16:13 UTC (rev 3441)
@@ -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.hh"
+
+///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
+
+}
+
Added: pkg/RcppSMC/src/simfunctions.cc
===================================================================
--- pkg/RcppSMC/src/simfunctions.cc (rev 0)
+++ pkg/RcppSMC/src/simfunctions.cc 2012-01-15 23:16:13 UTC (rev 3441)
@@ -0,0 +1,166 @@
+#include <iostream>
+#include <cmath>
+#include <gsl/gsl_randist.h>
+
+#include "smctc.hh"
+#include "simfunctions.hh"
+
+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));
+
+ while(y) {
+ lp += log(gsl_ran_ugaussian_pdf(y->value - x->value));
+ 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;
+}
+
Added: pkg/RcppSMC/src/simfunctions.hh
===================================================================
--- pkg/RcppSMC/src/simfunctions.hh (rev 0)
+++ pkg/RcppSMC/src/simfunctions.hh 2012-01-15 23:16:13 UTC (rev 3441)
@@ -0,0 +1,32 @@
+#include "smctc.hh"
+
+#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
More information about the Rcpp-commits
mailing list