geomc {geommc}

Markov chain Monte Carlo for discrete and continuous distributions using geometric MH algorithms.

Description

geomc produces Markov chain samples from a target distribution. The target can be a pdf or pmf. Users specify the target distribution by an R function that evaluates the log un-normalized pdf or pmf. geomc uses the geometric approach of Roy (2024) to move an uninformed base density (e.g. a random walk proposal) towards different global/local approximations of the target density. The base density can be specified along with its mean, covariance matrix, and a function for sampling from it. Gaussian densities can be specified by mean and variance only, although it is preferred to supply its density and sampling functions as well. If either or both of the mean and variance arguments and any of the density and sampling functions is missing, then a base density corresponding to a random walk with an appropriate scale parameter is used. One or more approximate target densities can be specified along with their means, covariance matrices, and a function for sampling from the densities. Gaussian densities can be specified by mean and variance only, although it is preferred to supply their densities and sampling functions as well. If either or both of the mean and variance arguments and any of the density and sampling functions is missing for the approximate target density, then a normal distribution with mean computed from a pilot run of a random walk Markov chain and a diagonal covariance matrix with a large variance is used. If the Argument gaus is set as FALSE then both the base and the approximate target can be specified by their densities and functions for sampling from it. That is, if gaus=FALSE, the functions specifying the means and variances of both the base and the approximate target densities are not used. If the target is a pmf (discrete distribution), then gaus=FALSE and imp [1]=TRUE (not the default values) need to be specified.

Usage

geomc(
  log.target,
  initial,
  n.iter,
  eps = 0.5,
  ind = FALSE,
  gaus = TRUE,
  imp = c(FALSE, n.samp = 1000, samp.base = FALSE),
  a = 1,
  mean.base,
  var.base,
  dens.base,
  samp.base,
  mean.ap.tar,
  var.ap.tar,
  dens.ap.tar,
  samp.ap.tar
)

Arguments

log.target

is the logarithm of the (un-normalized) target density function, needs to be written as a function of the current value x.

initial

is the initial values.

n.iter

is the no. of samples needed.

eps

is the value for epsilon perturbation. Default is 0.5.

ind

is False if either the base density, f or the approximate target density, g depends on the current value x. Default is False.

gaus

is True if both f and g are normal distributions. Default is True.

imp

is a vector of three elements. If gaus is TRUE, then the imp argument is not used. imp [1] is False if numerical integration is used, otherwise, importance sampling is used to compute \langle \sqrt{f}, \sqrt{g} \rangle. Default is False. imp [2] (n.samp) is no of samples in importance sampling. imp [3] (samp.base) is True if samples from f is used, otherwise samples from g is used. Default is False.

a

is the probability vector for the mixture proposal density. Default is the uniform distribution.

mean.base

is the mean of the base density f, needs to be written as a function of the current value x.

var.base

is the covariance matrix of the base density f, needs to be written as a function of the current value x.

dens.base

is the density function of the base density f, needs to be written as a function (y,x) (in this order) of the proposed value y and the current value x, although it may not depend on x.

samp.base

is the function to draw from the base density f, needs to be written as a function of the current value x.

mean.ap.tar

is the vector of means of the densities g_i(y|x), i=1,\dots,k. It needs to be written as a function of the current value x. It must have the same dimension as k times the length of initial.

var.ap.tar

is the matrix of covariance matrices of the densities g_i(y|x), i=1,\dots,k formed by column concatenation. It needs to be written as a function of the current value x. It must have the same dimension as the length of initial by k times the length of initial.

dens.ap.tar

is the vector of densities g_i(y|x), i=1,\dots,k. It needs to be written as a function (y,x) (in this order) of the proposed value y and the current value x, although it may not depend on x.

samp.ap.tar

is the function to draw from the densities g_i(y|x), i=1,\dots,k. It needs to be written as a function of (current value x, the indicator of mixing component kk). It must return a vector of the length of that of the initial.

Details

Using a geometric Metropolis-Hastings algorithm geom.mc produces Markov chains with the target as its stationary distribution. The details of the method can be found in Roy (2024).

Value

The function returns a list with the following elements:

samples

A matrix containing the MCMC samples. Each column is one sample.

acceptance.rate

The acceptance rate.

Author(s)

Vivekananda Roy [e-mail]

References

Roy, V.(2024) A geometric approach to informative MCMC sampling https://arxiv.org/abs/2406.09010

Examples

result <- geomc(log.target=function(y) dnorm(y,log=TRUE),initial=0,n.iter=500) 
#target is univariate normal
result$samples # the MCMC samples.
result$acceptance.rate # the acceptance rate.
result<-geomc(log.target=function(y) log(0.5*dnorm(y)+0.5*dnorm(y,mean=10,sd=1.4)),
initial=0,n.iter=500) #target is mixture of univariate normals, default choices
hist(result$samples)
result<-geomc(log.target=function(y) log(0.5*dnorm(y)+0.5*dnorm(y,mean=10,sd=1.4)),
initial=0,n.iter=500, mean.base = function(x) x,var.base= function(x) 4,
dens.base=function(y,x) dnorm(y, mean=x,sd=2),samp.base=function(x) x+2*rnorm(1),
mean.ap.tar=function(x) c(0,10),var.ap.tar=function(x) c(1,1.4^2),
dens.ap.tar=function(y,x) c(dnorm(y),dnorm(y,mean=10,sd=1.4)),
samp.ap.tar=function(x,kk=1){if(kk==1){return(rnorm(1))} else{return(10+1.4*rnorm(1))}})
#target is mixture of univariate normals, random walk base density, an informed 
#choice for dens.ap.tar
hist(result$samples)
samp.ap.tar=function(x,kk=1){s.g=sample.int(2,1,prob=c(.5,.5))
if(s.g==1){return(rnorm(1))
}else{return(10+1.4*rnorm(1))}}
result<-geomc(log.target=function(y) log(0.5*dnorm(y)+0.5*dnorm(y,mean=10,sd=1.4)),
initial=0,n.iter=500,gaus=FALSE,imp=c(TRUE,n.samp=100,samp.base=TRUE),
dens.base=function(y,x) dnorm(y, mean=x,sd=2),samp.base=function(x) x+2*rnorm(1),
dens.ap.tar=function(y,x) 0.5*dnorm(y)+0.5*dnorm(y,mean=10,sd=1.4),
samp.ap.tar=samp.ap.tar)
#target is mixture of univariate normals, random walk base density, another 
#informed choice for dens.ap.tar
hist(result$samples)
result <- geomc(log.target=function(y) -0.5*crossprod(y),initial=rep(0,4),
n.iter=500) #target is multivariate normal, default choices
rowMeans(result$samples)
size=5
result <- geomc(log.target = function(y) dbinom(y, size, 0.3, log = TRUE),
initial=0,n.iter=500,ind=TRUE,gaus=FALSE,imp=c(TRUE,n.samp=1000,samp.base=TRUE),
dens.base=function(y,x) 1/(size+1), samp.base= function(x) sample(seq(0,size,1),1),
dens.ap.tar=function(y,x) dbinom(y, size, 0.7),samp.ap.tar=function(x,kk=1) rbinom(1, size, 0.7))
 #target is binomial
 table(result$samples)