SAMCPLUS: SAMC Sampler with C++

In SAMCpack: Stochastic Approximation Monte Carlo (SAMC) Sampler and Methods

Description

The function SAMCPLUS is a generic SAMC sampler for distributions on continuous domain. Instead of an R function, SAMCPLUS requires a function pointer to be provided for faster sampling, with all other values and parameters being equal to its cousin SAMC. We limited the flexibility of the function pointer to be passed. See the below for more details or the vignette.

Usage

SAMCPLUS(nv, energy, data = NA, options = list())

Arguments

nv	number of variables.
energy	a CPP function pointer for negative log density.
data	extra data to be supplemented. It should be a vector, a matrix, or a list.
options	a list specifying parameters/options for SAMC algorithm, PARAMETER SPECIFICATION DESCRIPTION domain vector(2) or matrix((nv\times 2)) domain of sample space partition vector(m) energy partition vecpi vector(m-1) desired sampling distribution tau positive number temperature niter positive integer number of iterations to be run t0 (0.5,1] gain factor sequence value xi positive number gain factor sequence exponent stepsize positive number stepsize for random-walk sampler trange vector(2) or matrix(m\times 2) domain of estimates for \log(g_i /π_i)

Value

a named list containing

samples

an (niter\times nv) samples generated.

frequency

length-m vector of visiting frequency for energy partition.

theta

length-m vector of estimates of \log(g_i / π_i)

Note on writing your own C++ function

First, the output should be returned as SEXP rather than double in evaluating the negative log density. Second, the variable and extra data should be provided as arma::vec and arma::mat type, with an exception for Rcpp::List for list-valued data. This means, for the data Even though we could let data to be passed freely, we believe using RcppArmadillo, which is a templated C++ linear algebra library, enables easier writing of one's own C++ code in a style of R or MATLAB while providing sufficient computational power. Furthermore, limiting extra data to one of 3 types (vector, matrix, and list) reduces potential type-matching issue in encapsulating of the current environment by removing unexpected errors a user might have incurred.

Author(s)

Kisung You

References

\insertRef

SAMCSAMCpack

Examples

##### Two-Dimensional Multimodal sampling
## Step 1 : Define negative log-density function as a CPP function
cppscript = "SEXP funcH(arma::vec x){
double x1 = x(0);
double x2 = x(1);
double val1 = (-std::pow((x1*sin(20*x2)+x2*sin(20*x1)),2))*cosh(sin(10*x1)*x1);
double val2 = (-std::pow((x1*cos(10*x2)-x2*sin(10*x1)),2))*cosh(cos(20*x2)*x2);
return Rcpp::wrap(val1+val2);
}"
func_ptr = RcppXPtrUtils::cppXPtr(cppscript,depends="RcppArmadillo") # as a pointer

## Step 2 : Prepare a setting option
myoption = list()
myoption$partition = c(-Inf,seq(from=-8,to=0,length.out=41))
myoption$tau       = 1.0
myoption$domain    = c(-1.1,1.1)
myoption$vecpi     = as.vector(rep(1/41,41))
myoption$niter     = 200001
myoption$stepsize  = 0.25

## Step 3 : Run The Code
res = SAMCPLUS(2,func_ptr,options=myoption)

## Step 4 : Visualize
select = seq(from=101,to=myoption$niter,by=100) # 100 burn-in, 1/100 thinning 
par(mfrow=c(1,2))
plot(res$samples[select,1], res$samples[select,2],xlab='x',ylab='y',main='samples')
barplot(as.vector(res$frequency/sum(res$frequency)),
        main="visiting frequency by energy partition",
        names.arg=myoption$partition[-1], xlab="energy")

##### (2) Use Extra Data
## Define negative log-density function as CPP function
cppscript2 = "SEXP funcH2(arma::vec x, arma::vec data){
double x1 = x(0);
double x2 = x(1);
double par1 = data(0);
double par2 = data(1);

double val1 = (-std::pow((x1*sin(par2*x2)+x2*sin(par2*x1)),2))*cosh(sin(par1*x1)*x1);
double val2 = (-std::pow((x1*cos(par1*x2)-x2*sin(par1*x1)),2))*cosh(cos(par2*x2)*x2);
return Rcpp::wrap(val1+val2);
}"
func_ptr2 = RcppXPtrUtils::cppXPtr(cppscript2,depends="RcppArmadillo") # as a pointer

## Run The Code
vecdata = as.vector(c(10,20)) 
res2 = SAMCPLUS(2,func_ptr2,data=vecdata, options=myoption)
select = seq(from=101,to=ex_niter,by=100) # 100 burn-in, 1/100 thinning 
par(mfrow=c(1,2))
plot(res2$samples[select,1], res2$samples[select,2],xlab='x',ylab='y',main='samples')
barplot(as.vector(res2$frequency/sum(res2$frequency)),
        main="visiting frequency by energy partition",
        names.arg=ex_part[2:(m+1)], xlab="energy")

SAMCpack documentation built on May 2, 2019, 7:31 a.m.