roll minor version and date

also push a minor update to an old example made a few days ago

Author: eddelbuettel

ChangeLog

@@ -1,3 +1,8 @@
+2020-04-01  Dirk Eddelbuettel  <[email protected]>
+
+	* DESCRIPTION (Version, Date): Roll minor version
+        * inst/include/Rcpp/config.h: Idem
+
 2020-03-31  Dirk Eddelbuettel  <[email protected]>
 
 	* docker/ci-dev/Dockerfile: Also install the 'codetools' package
@@ -6,6 +11,10 @@
 
 	* R/Attributes.R: Fix for sourceCpp() on Windows w/R-devel
 
+2020-03-28  Dirk Eddelbuettel  <[email protected]>
+
+	* inst/examples/RcppGibbs/RcppGibbs_Updated.R: Updated example
+
 2020-03-24  Dirk Eddelbuettel  <[email protected]>
 
 	* DESCRIPTION (Version, Date): Roll minor version

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: Rcpp
 Title: Seamless R and C++ Integration
-Version: 1.0.4.5
-Date: 2020-03-24
+Version: 1.0.4.6
+Date: 2020-04-01
 Author: Dirk Eddelbuettel, Romain Francois, JJ Allaire, Kevin Ushey, Qiang Kou,
  Nathan Russell, Douglas Bates and John Chambers
 Maintainer: Dirk Eddelbuettel <[email protected]>

inst/examples/RcppGibbs/RcppGibbs_Updated.R

@@ -0,0 +1,189 @@
+## Simple Gibbs Sampler Example
+## Adapted from Darren Wilkinson's post at
+## http://darrenjw.wordpress.com/2010/04/28/mcmc-programming-in-r-python-java-and-c/
+##
+## Sanjog Misra and Dirk Eddelbuettel, June-July 2011
+## Updated by Dirk Eddelbuettel, Mar 2020
+
+suppressMessages({
+    library(Rcpp)
+    library(rbenchmark)
+})
+
+
+## Actual joint density -- the code which follow implements
+## a Gibbs sampler to draw from the following joint density f(x,y)
+fun <- function(x,y) {
+    x*x * exp(-x*y*y - y*y + 2*y - 4*x)
+}
+
+## Note that the full conditionals are propotional to
+## f(x|y) = (x^2)*exp(-x*(4+y*y))              : a Gamma density kernel
+## f(y|x) = exp(-0.5*2*(x+1)*(y^2 - 2*y/(x+1)) : Normal Kernel
+
+## There is a small typo in Darrens code.
+## The full conditional for the normal has the wrong variance
+## It should be 1/sqrt(2*(x+1)) not 1/sqrt(1+x)
+## This we can verify ...
+## The actual conditional (say for x=3) can be computed as follows
+## First - Construct the Unnormalized Conditional
+fy.unnorm <- function(y) fun(3,y)
+
+## Then - Find the appropriate Normalizing Constant
+K <- integrate(fy.unnorm,-Inf,Inf)
+
+## Finally - Construct Actual Conditional
+fy <- function(y) fy.unnorm(y)/K$val
+
+## Now - The corresponding Normal should be
+fy.dnorm <- function(y) {
+    x <- 3
+    dnorm(y,1/(1+x),sqrt(1/(2*(1+x))))
+}
+
+## and not ...
+fy.dnorm.wrong <- function(y) {
+    x <- 3
+    dnorm(y,1/(1+x),sqrt(1/((1+x))))
+}
+
+if (interactive()) {
+    ## Graphical check
+    ## Actual (gray thick line)
+    curve(fy,-2,2,col='grey',lwd=5)
+
+    ## Correct Normal conditional (blue dotted line)
+    curve(fy.dnorm,-2,2,col='blue',add=T,lty=3)
+
+    ## Wrong Normal (Red line)
+    curve(fy.dnorm.wrong,-2,2,col='red',add=T)
+}
+
+## Here is the actual Gibbs Sampler
+## This is Darren Wilkinsons R code (with the corrected variance)
+## But we are returning only his columns 2 and 3 as the 1:N sequence
+## is never used below
+Rgibbs <- function(N,thin) {
+    mat <- matrix(0,ncol=2,nrow=N)
+    x <- 0
+    y <- 0
+    for (i in 1:N) {
+        for (j in 1:thin) {
+            x <- rgamma(1,3,y*y+4)
+            y <- rnorm(1,1/(x+1),1/sqrt(2*(x+1)))
+        }
+        mat[i,] <- c(x,y)
+    }
+    mat
+}
+
+## Now for the Rcpp version -- Notice how easy it is to code up!
+
+cppFunction("NumericMatrix RcppGibbs(int N, int thn){
+    NumericMatrix mat(N, 2);        // Setup storage
+    double x = 0, y = 0;            // The rest follows the R version
+    for (int i = 0; i < N; i++) {
+        for (int j = 0; j < thn; j++) {
+            x = R::rgamma(3.0,1.0/(y*y+4));
+            y = R::rnorm(1.0/(x+1),1.0/sqrt(2*x+2));
+        }
+        mat(i,0) = x;
+        mat(i,1) = y;
+    }
+    return mat;             // Return to R
+}")
+
+
+## Use of the sourceCpp() is preferred for users who wish to source external
+## files or specify their headers and Rcpp attributes within their code.
+## Code here is able to easily be extracted and placed into its own C++ file.
+
+## Compile and Load
+sourceCpp(code="
+#include <RcppGSL.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+using namespace Rcpp;  // just to be explicit
+
+// [[Rcpp::depends(RcppGSL)]]
+
+// [[Rcpp::export]]
+NumericMatrix GSLGibbs(int N, int thin){
+    gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+    double x = 0, y = 0;
+    NumericMatrix mat(N, 2);
+    for (int i = 0; i < N; i++) {
+        for (int j = 0; j < thin; j++) {
+            x = gsl_ran_gamma(r,3.0,1.0/(y*y+4));
+            y = 1.0/(x+1)+gsl_ran_gaussian(r,1.0/sqrt(2*x+2));
+        }
+        mat(i,0) = x;
+        mat(i,1) = y;
+    }
+    gsl_rng_free(r);
+
+    return mat;           // Return to R
+}")
+
+
+
+## Now for some tests
+## You can try other values if you like
+## Note that the total number of interations are N*thin!
+Ns <- c(1000,5000,10000,20000)
+thins <- c(10,50,100,200)
+tim_R <- rep(0,4)
+tim_Rgsl <- rep(0,4)
+tim_Rcpp <- rep(0,4)
+
+for (i in seq_along(Ns)) {
+    tim_R[i] <- system.time(mat <- Rgibbs(Ns[i],thins[i]))[3]
+    tim_Rgsl[i] <- system.time(gslmat <- GSLGibbs(Ns[i],thins[i]))[3]
+    tim_Rcpp[i] <- system.time(rcppmat <- RcppGibbs(Ns[i],thins[i]))[3]
+    cat("Replication #", i, "complete \n")
+}
+
+## Comparison
+speedup <- round(tim_R/tim_Rcpp,2);
+speedup2 <- round(tim_R/tim_Rgsl,2);
+summtab <- round(rbind(tim_R, tim_Rcpp,tim_Rgsl,speedup,speedup2),3)
+colnames(summtab) <- c("N=1000","N=5000","N=10000","N=20000")
+rownames(summtab) <- c("Elasped Time (R)","Elapsed Time (Rcpp)", "Elapsed Time (Rgsl)",
+                       "SpeedUp Rcpp", "SpeedUp GSL")
+print(summtab)
+
+## Contour Plots -- based on Darren's example
+if (interactive() && require(KernSmooth)) {
+    op <- par(mfrow=c(4,1),mar=c(3,3,3,1))
+    x <- seq(0,4,0.01)
+    y <- seq(-2,4,0.01)
+    z <- outer(x,y,fun)
+    contour(x,y,z,main="Contours of actual distribution",xlim=c(0,2), ylim=c(-2,4))
+    fit <- bkde2D(as.matrix(mat),c(0.1,0.1))
+    contour(drawlabels=T, fit$x1, fit$x2, fit$fhat, xlim=c(0,2), ylim=c(-2,4),
+            main=paste("Contours of empirical distribution:",round(tim_R[4],2)," seconds"))
+    fitc <- bkde2D(as.matrix(rcppmat),c(0.1,0.1))
+    contour(fitc$x1,fitc$x2,fitc$fhat,xlim=c(0,2), ylim=c(-2,4),
+            main=paste("Contours of Rcpp based empirical distribution:",round(tim_Rcpp[4],2)," seconds"))
+    fitg <- bkde2D(as.matrix(gslmat),c(0.1,0.1))
+    contour(fitg$x1,fitg$x2,fitg$fhat,xlim=c(0,2), ylim=c(-2,4),
+            main=paste("Contours of GSL based empirical distribution:",round(tim_Rgsl[4],2)," seconds"))
+    par(op)
+}
+
+
+## also use rbenchmark package
+N <- 20000
+thn <- 200
+res <- benchmark(Rgibbs(N, thn),
+                 RcppGibbs(N, thn),
+                 GSLGibbs(N, thn),
+                 columns=c("test", "replications", "elapsed",
+                           "relative", "user.self", "sys.self"),
+                 order="relative",
+                 replications=10)
+print(res)
+
+
+## And we are done

inst/include/Rcpp/config.h

@@ -30,7 +30,7 @@
 #define RCPP_VERSION_STRING     "1.0.4"
 
 // the current source snapshot
-#define RCPP_DEV_VERSION        RcppDevVersion(1,0,4,5)
-#define RCPP_DEV_VERSION_STRING "1.0.4.5"
+#define RCPP_DEV_VERSION        RcppDevVersion(1,0,4,6)
+#define RCPP_DEV_VERSION_STRING "1.0.4.6"
 
 #endif