Particle marginal Metropolis Hastings Example

ChangeLog

@@ -1,3 +1,19 @@
+2017-08-03  Leah South  <leah.south@hdr.qut.edu.au>
+
+	* inst/include/sampler.h: Initialising dlogNCPath to 0.
+	* src/nonLinPMMH.cpp: A new example based on using PMMH for a
+	non-linear state space model
+	* inst/include/nonLinPMMH.h: Idem.
+	* R/nonLinPMMH.R: Idem.
+	* man/nonLinPMMH.Rd: Documenting the new example.
+	* NAMESPACE: Adding the new function nonLinPMMH.
+
+	* R/simNonlin.R: Generalising the function so that it can be
+	used for both nonLinPMMH and pfNonlinBS.
+	* man/simNonlin.Rd: Documenting the function separately because
+	it is used in multiple examples.
+	* man/pfNonlinBS.Rd: Documenting simNonlin separately.
+
 2017-08-02  Dirk Eddelbuettel  <edd@debian.org>
 
 	* R/LinReg.R: Silence R CMD check --as-cran about data(radiate)

NAMESPACE

@@ -4,7 +4,7 @@ useDynLib("RcppSMC", .registration=TRUE)
 import("methods")
 importFrom("Rcpp", "evalCpp")               # with Rcpp 0.11.0
 importFrom("grDevices", "dev.new")
-importFrom("graphics", "lines", "par", "plot", "points", "polygon", "title")
+importFrom("graphics", "lines", "par", "plot", "points", "polygon", "title", "abline")
 importFrom("stats", "rnorm", "rt")
 importFrom("utils", "read.table", "data")
 
@@ -14,6 +14,7 @@ export("blockpfGaussianOpt",
        "pfNonlinBS",
        "LinReg",
        "LinRegLA",
+       "nonLinPMMH",
        "simGaussian",
        "simLineart",
        "simNonlin")

R/RcppExports.R

@@ -13,6 +13,10 @@ LinRegLA_impl <- function(Data, intemps, lNumber) {
     .Call('RcppSMC_LinRegLA_impl', PACKAGE = 'RcppSMC', Data, intemps, lNumber)
 }
 
+nonLinPMMH_impl <- function(data, lNumber, lMCMCits) {
+    .Call('RcppSMC_nonLinPMMH_impl', PACKAGE = 'RcppSMC', data, lNumber, lMCMCits)
+}
+
 pfLineartBS_impl <- function(data, part, usef, fun) {
     .Call('RcppSMC_pfLineartBS_impl', PACKAGE = 'RcppSMC', data, part, usef, fun)
 }

R/nonLinPMMH.R

@@ -0,0 +1,45 @@
+nonLinPMMH<- function(data, particles=5000, iterations=10000, burnin = 0, plot=FALSE) {
+
+    if (missing(data)) {
+         warning("data argument contained no data, using data simulated from the model.")
+         data <- simNonlin(len=500,var_init=5,var_evol=10,var_obs=1,cosSeqOffset=0)$data
+    }
+
+    if (burnin>=iterations) {
+        stop("Burn-in must be less than iterations.")
+    }
+
+    res <- nonLinPMMH_impl(as.matrix(data), particles, iterations)
+
+	res.plot <- res[burnin+1:iterations,]
+
+    # Replicating Figure 4(c) from Andrieu et al. (2010).
+    # May want to add autocorrelation plots of the parameters, like in Figures 5(b) and 5(d).
+    if (plot) {
+        par(mfrow=c(2,3),oma=c(0,0,2,0))
+        with(res.plot, hist(samples_sigv,
+                       xlab=expression(sigma_v),ylab='density',
+                       main = NA))
+        abline(v=sqrt(10),lty=2, col='darkblue',lwd=2)
+        with(res.plot, plot(samples_sigw,samples_sigv,
+                       xlab=expression(sigma_w),ylab=expression(sigma_v),
+                       main = NA, col='darkblue'))
+        with(res.plot, plot(samples_sigv,type="l",
+                       ylab=expression(sigma_v),
+                       main = NA,xlim=c(0,iterations-burnin)))
+        with(res.plot, plot(samples_sigv,samples_sigw,
+                       xlab=expression(sigma_v),ylab=expression(sigma_w),
+                       main = NA, col='darkblue'))
+        with(res.plot, hist(samples_sigw,
+                       xlab=expression(sigma_w),ylab='density',
+                       main = NA))
+        abline(v=1,lty=2, col='darkblue',lwd=2)
+        with(res.plot, plot(samples_sigw,type="l",
+                       ylab=expression(sigma_w),
+                       main = NA,xlim=c(0,iterations-burnin)))
+        title("Posterior Estimates",outer=TRUE)
+        par(mfrow=c(1,1))
+    }  
+
+    invisible(res)
+}

R/simNonlin.R

@@ -1,14 +1,14 @@
-simNonlin <- function(len = 50)
+simNonlin <- function(len = 50, var_init = 10, var_evol = 10, var_obs = 1, cosSeqOffset = -1)
 {
    sim <- list()
 
-   innovations <- rnorm(len) * sqrt(10)
-   sim$state[1] <- innovations[1]
+   sim$state[1] <- rnorm(1) * sqrt(var_init)
+   innovations <- rnorm(len-1) * sqrt(var_evol)
    for (i in 2:len) {
        sim$state[i] <- 0.5 * sim$state[i-1] + 25 * sim$state[i-1] /
-           (1 + sim$state[i-1]^2) + 8 * cos(1.2*(i-1)) + innovations[i]
+           (1 + sim$state[i-1]^2) + 8 * cos(1.2*(i+cosSeqOffset)) + innovations[i-1]
    }
-   sim$data <- sim$state^2 / 20 + rnorm(len)
+   sim$data <- sim$state^2 / 20 + rnorm(len) * sqrt(var_obs)
 
    invisible(sim)
 }

inst/include/nonLinPMMH.h

@@ -0,0 +1,41 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
+//
+// nonLinPMMH.h: Example 3.1 of Andrieu et al. (2010). Implementing particle marginal
+// Metropolis-Hastings for a toy non-linear state space model previously described in
+// Gordon et al. (1993) and Kitagawa (1996).
+//
+// Copyright (C) 2017         Dirk Eddelbuettel, Adam Johansen and Leah South
+//
+// 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"
+
+namespace nonLinPMMH {
+
+    class parameters 
+    {
+    public:
+        double sigv, sigw;
+    };
+    arma::vec y; //data
+
+    double logPrior(const parameters & proposal);
+    void fInitialise(double & value, double & logweight);
+    void fMove(long lTime, double & value, double & logweight);
+
+    parameters theta_prop;
+}
+

inst/include/sampler.h

@@ -232,6 +232,7 @@ namespace smc {
     void sampler<Space>::Initialise(void)
     {
         T = 0;
+        dlogNCPath = 0.0;
 
         //Set the initial values and log weights of the particles
         std::vector<Space> InitVal(N);

man/nonLinPMMH.Rd

@@ -0,0 +1,74 @@
+\name{nonLinPMMH}
+\alias{nonLinPMMH}
+\title{Particle marginal Metropolis-Hastings for a non-linear state space model.}
+\usage{
+nonLinPMMH(data, particles = 5000, iterations = 10000, burnin = 0,
+  plot = FALSE)
+}
+\arguments{
+\item{data}{A vector of the observed data.}
+
+\item{particles}{An integer specifying the number of particles in the particle
+filtering estimates of the likelihood.}
+
+\item{iterations}{An integer specifying the number of MCMC iterations.}
+
+\item{burnin}{The number of iterations to remove from the beginning of the MCMC chain
+(for plotting purposes only).}
+
+\item{plot}{A boolean variable to determine whether to plot the posterior estimates
+and MCMC chain.}
+}
+\value{
+A \code{data.frame} containing the chain
+of simulated \eqn{\sigma_v} and \eqn{\sigma_w} values, as well as the
+corresponding log likelihood estimates and log prior values.
+}
+\description{
+The \code{nonLinPMMH} function implements particle marginal Metropolis Hastings
+for the non-linear state space model described in Section 3.1 of
+Andrieu et al. (2010).
+}
+\details{
+This example uses particle marginal Metropolis Hastings to estimate
+the standard deviation of the evolution and observation noise in the following
+non-linear state space model:
+
+	\eqn{x(n) = 0.5 x(n-1) + 25 x(n-1) / (1+x(n-1)^2) + 8 cos(1.2n)+ e(n)} and 
+
+	\eqn{y(n) = x(n)^2 / 20 + f(n)}
+
+where e(n) and f(n) are mutually-independent normal random
+variables of variances var_evol and var_obs, respectively,
+and \eqn{x(0) ~ N(0,5)}.
+
+Following Andrieu, Doucet and Holenstein (2010), the priors are
+\eqn{var_evol ~ IG(0.01,0.01)} and \eqn{var_obs ~ IG(0.01,0.01)} where IG
+is the inverse gamma distribution.
+
+Data can be simulated from the model using \code{\link{simNonlin}}.
+}
+\examples{
+\dontrun{
+sim <- simNonlin(len=500,var_init=5,var_evol=10,var_obs=1,cosSeqOffset=0)
+res <- nonLinPMMH(sim$data,particles=5000,iterations=50000,burnin=10000,plot=TRUE)
+}
+
+}
+\references{
+C. Andrieu, A. Doucet, and R. Holenstein. Particle Markov chain Monte Carlo methods.
+Journal of the Royal Statistical Society: Series B (Statistical Methodology),
+72(3):269-342, 2010.
+}
+\seealso{
+\code{\link{simNonlin}} for a function to simulate from the model and
+\code{\link{pfNonlinBS}} for a simple bootrap particle filter
+applied to a similar non-linear state space model.
+}
+\author{
+Adam M. Johansen, Dirk Eddelbuettel and Leah F. South
+}
+\concept{
+Bayesian PMMH PF
+}
+\keyword{programming}

man/pfNonlinBS.Rd

@@ -1,32 +1,25 @@
 \name{pfNonlinBS}
 \alias{pfNonlinBS}
-\alias{simNonlin}
 \title{Nonlinear Bootstrap Particle Filter (Univariate Non-Linear State Space Model)}
 \description{
   The \code{pfNonlinBS} function provides a simple example for
   \pkg{RcppSMC}. It is a simple \dQuote{bootstrap} particle filter which employs
   multinomial resampling after each iteration applied to the ubiquitous "nonlinear
   state space model" following Gordon, Salmond and Smith (1993).
-
-  The \code{simNonlin} function simulates data from the associated model.
 }
 \usage{
   pfNonlinBS(data, particles=500, plot=FALSE) 
-  simNonlin(len)
 }
 \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.} 
-  \item{len}{The length of data sequence to simulate.}
+    standard deviation intervals.}
 }
 \value{
   The \code{pfNonlinBS} function returns two vectors, the first containing the posterior
   filtering means; the second the posterior filtering standard deviations.
-
-  The \code{simNonlin} function returns a list containing the state and data sequences. 
 }
 \details{
   The \code{pfNonlinbs} function provides a simple example for
@@ -37,9 +30,7 @@
   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. 
-
-  The \code{simNonlin} function simulates from the same model.
+  resampling after each iteration.
 }
 \references{
   N. J. Gordon, S. J. Salmond, and A. F. M. Smith. Novel approach to
@@ -50,5 +41,10 @@
   sim <- simNonlin(len=50)
   res <- pfNonlinBS(sim$data,particles=500,plot=TRUE)
 }
-\author{Adam M. Johansen and Dirk Eddelbuettel}
+\seealso{
+\code{\link{simNonlin}} for a function to simulate from the model and
+\code{\link{nonLinPMMH}} for an example of particle
+marginal Metropolis Hastings applied to a non-linear state space model.
+}
+\author{Adam M. Johansen, Dirk Eddelbuettel and Leah F. South}
 \keyword{programming}

man/simNonlin.Rd

@@ -0,0 +1,66 @@
+\name{simNonlin}
+\alias{simNonlin}
+\title{Simulates from a simple nonlinear state space model.}
+\usage{
+simNonlin(len = 50, var_init = 10, var_evol = 10, var_obs = 1,
+  cosSeqOffset = -1)
+}
+\arguments{
+\item{len}{The length of data sequence to simulate.}
+
+\item{var_init}{The variance of the noise for the initial state.}
+
+\item{var_evol}{The variance of the noise for the state evolution .}
+
+\item{var_obs}{The variance of the observation noise.}
+
+\item{cosSeqOffset}{This is related to the indexing in the
+cosine function in the evoluation equation. A value of -1
+can be used to follow the specification of Gordon, Salmond
+and Smith (1993) and 0 can be used to follow
+Andrieu, Doucet and Holenstein (2010).}
+}
+\value{
+The \code{simNonlin} function returns a list containing the state and data sequences.
+}
+\description{
+The \code{simNonlin} function simulates data from the models used
+in \code{link{pfNonlinBS}} and \code{link{nonLinPMMH}}.
+}
+\details{
+The \code{simNonlin} function simulates from
+a simple nonlinear state space model with
+state evolution and observation equations:
+
+	\eqn{x(n) = 0.5 x(n-1) + 25 x(n-1) / (1+x(n-1)^2) + 8 cos(1.2(n+cosSeqOffset))+ e(n)} and 
+
+	\eqn{y(n) = x(n)^2 / 20 + f(n)}
+
+where \eqn{e(n)} and \eqn{f(n)} are mutually-independent normal random
+variables of variances var_evol and var_obs, respectively,
+and \eqn{x(0) ~ N(0,var_init)}.
+
+Different variations of this model can be found in
+Gordon, Salmond and Smith (1993) and
+Andrieu, Doucet and Holenstein (2010). A cosSeqOffset
+of -1 is consistent with the former and
+0 is consistent with the latter.
+}
+\references{
+C. Andrieu, A. Doucet, and R. Holenstein. Particle Markov chain Monte Carlo methods.
+Journal of the Royal Statistical Society: Series B (Statistical Methodology),
+72(3):269-342, 2010.
+
+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.
+}
+\author{
+Adam M. Johansen, Dirk Eddelbuettel and Leah F. South
+}
+\seealso{
+\code{\link{pfNonlinBS}} for a simple bootrap particle filter
+applied to this model and \code{\link{nonLinPMMH}} for particle
+marginal Metropolis Hastings applied to estimating the standard
+deviation of the state evolution and observation noise.
+}

src/RcppExports.cpp

@@ -44,6 +44,19 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// nonLinPMMH_impl
+Rcpp::DataFrame nonLinPMMH_impl(arma::vec data, unsigned long lNumber, unsigned long lMCMCits);
+RcppExport SEXP RcppSMC_nonLinPMMH_impl(SEXP dataSEXP, SEXP lNumberSEXP, SEXP lMCMCitsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::vec >::type data(dataSEXP);
+    Rcpp::traits::input_parameter< unsigned long >::type lNumber(lNumberSEXP);
+    Rcpp::traits::input_parameter< unsigned long >::type lMCMCits(lMCMCitsSEXP);
+    rcpp_result_gen = Rcpp::wrap(nonLinPMMH_impl(data, lNumber, lMCMCits));
+    return rcpp_result_gen;
+END_RCPP
+}
 // pfLineartBS_impl
 Rcpp::DataFrame pfLineartBS_impl(arma::mat data, unsigned long part, bool usef, Rcpp::Function fun);
 RcppExport SEXP RcppSMC_pfLineartBS_impl(SEXP dataSEXP, SEXP partSEXP, SEXP usefSEXP, SEXP funSEXP) {