SAMCpack (version 0.1.1)

SAMCPLUS: SAMC Sampler with C++

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

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 / \pi_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.

References

SAMCSAMCpack

Examples

# NOT RUN {
##### 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")
# }
# NOT RUN {
# }