SAMC: Stochastic Approximation Monte Carlo (SAMC) Sampler

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

Description

The function SAMC is a generic SAMC sampler for distributions on continuous domain. An R function for negative log density of your choice is required, as well as some parameters in SAMC framework.

Usage

SAMC(nv, energy, options = list())

Arguments

nv	number of variables.
energy	an R function for negative log density.
options	a list specifying parameters/options for SAMC algorithm. Below, vector(k) means a vector of length k, and matrix likewise. 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)

Author(s)

Kisung You

References

\insertRef

SAMCSAMCpack

Examples

##### Two-Dimensional Multimodal sampling
## Step 1 : Define negative log-density function as an R function
func_r = function(x){
x1 = x[1]; x2 = x[2];
val1 = (-(x1*sin(20*x2)+x2*sin(20*x1))^2)*cosh(sin(10*x1)*x1);
val2 = (-(x1*cos(10*x2)-x2*sin(10*x1))^2)*cosh(cos(20*x2)*x2);
return(val1+val2);
}

## 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     = 20000
myoption$stepsize  = c(0.25, 10)

## Step 3 : Run The Code
res = SAMC(2,func_r,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")

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