R: Metropolis Hasting Algorithm Constrained on a Simplex

RunMh {SALTSampler}

R Documentation

Metropolis Hasting Algorithm Constrained on a Simplex

Description

This function runs the Metropolis Hasting algorithm constrained on a simplex. The function can be used with any target distribution on the simplex defined by the user. Alternatively, two common target distributions are built into the function and can be specifed by the user. The function is designed to continue to perform well in difficult cases, such as those in high dimensions or with parameters that differ by orders of magnitude. Care is also taken to ensure accuracy even when some coordinates are numerically close to 0 or 1.

Usage

RunMh(center, B, concentration = 1, h, type = 'user', dat = NULL, pars = NULL)

Arguments

`center`	Vector of numeric values summing to 1 that define the center of the distributional parameters of the posterior. For type `'dirichlet'`, the parameter `a` is defined such that `a_{i}` is the ith element of `center` times `concentration`. For type `'multinom'`, the multinomial distribution parameter, `p_{i}`, is the ith value of `center`
`B`	Number of iterations to run the chain
`concentration`	This argument specifies the concentration parameter where `a` is defined such that `a_{i}` is the ith element of `center` times `concentration`. This is typically used with type `'dirichlet'`, but can also be used in a user-defined function. This arguments defaults to 1, so has no effect if it is not specified.
`h`	Vector of step sizes. Length of vector must match length of `center`
`type`	Specifies the target distribution. Select type `'user'` if a target distribution has already been defined (see details). Select type `'dirichlet'` for a Dirichlet distribution and type `'multinom'` for a multinomial distribution
`dat`	A matrix or vector passing data to the sampler. For type `'multinom'`, this is a matrix giving data from repeated multinomial draws where the data is formatted in the same way as data obtained via `GenData`. The number of the items in the ith bin on the jth multinomial trial should be in the ith column and the jth row of the matrix. For type `'user'`, any matrix or vector of data can be used to match the form specified in the user's target function. If unspecified, this argument defaults to `NULL`
`pars`	A list of additional parameters that can be passed to the user-specified target function for type `'user'` if desired. Argument defaults to `NULL`

Details

Any target distribution on the simplex can be used with this function by defining a target distribution function in the environment prior to running RunMh. The function should be named Target and should take in parameters ycand and ycurrent, which are the current and proposed samples on the logit scale, and parameter a, which is center times concentration. Parameters dat and pars can be set to NULL. Alternatively, dat can be used to provide data to the target function and/or pars can be used to provide a list of additional parameters to the the target function. The target function should output the ratio of the log-likelihood of the posterior distribution for the proposal, \theta = ycand, to the log-likelihood of the posterior for the current value, \theta = ycurrent. For simple cases, there are built-in target distributions. For type 'dirichlet', RunMh uses a Dirichlet distribution as a posterior distribution. For type 'multinomial', RunMh samples the distributional parameters of a multinomial distribution that would have generated the data inputted for dat.

Value

An object of class mhOut. mhOut has 12 attributes.

`Y`	Matrix of MCMC samples on logit scale
`S`	Matrix of MCMC samples on true scale
`runTime`	Summary of the MCMC runtime. The first entry gives the total user CPU time, the second entry gives the system CPU time, and the third entry gives the true elapsed time
`moveCount`	Number of steps where the proposal value was accepted
`p`	Length of `center` vector
`center`	Vector of numeric values summing to `1` that help to define distributional parameters. For type `'dirichlet'`, the parameter `a` is defined such that `a_{i}` is the ith element of `center` times `concentration`. For type `'multinom'`, the multinomial distribution parameter, `p_{i}`, is the ith value of `center`
`B`	Number of iterations to run the chain
`concentration`	For type `'dirichlet'`, this argument specifies the concentration parameter where `a` is defined such that `a_{i}` is the ith element of `center` times `concentration`. Otherwise, this argument takes on its default value of 1 and has no effect
`h`	Vector of step sizes. Length of vector must match length of `center`
`type`	Specifies the target distribution. Select type `'user'` if a target distribution has already been defined (see details). Select type `'dirichlet'` for a Dirichlet distribution and type `'multinom'` for a multinomial distribution
`dat`	A matrix or vector passing data to the sampler. For type `'multinom'`, a matrix giving data from repeated multinomial draws where the data is formatted in the same way as data obtained via `GenData`. The number of the items in the ith bin on the jth multinomial trial should be in the ith column and the jth row of the matrix. For type `'user'`, any matrix or vector of data can be used to match the form specified in the user's target function. If unspecified, this argument defaults to `NULL`
`a`	Dirichlet distribution parameters, `a`, where `a_{i}`, is the ith element of `center` times `concentration`

Examples

###Dirichlet sampling in 3-simplex
dir <- RunMh(center = c(0.7, 0.2, 0.1), B = 2e3, concentration = 10,
                        h = c(2, 2, 2), type = 'dirichlet', dat = NULL)
                        
####Multinomial sampling                  
## Not run: 
sampData <- GenData(center = c(0.2, 0.3, 0.5), n = 100, size = 10)
multinom <- RunMh(center = c(0.2, 0.3, 0.5), B = 1e4, h = c(2,2,2), 
                  type = 'multinom', dat = sampData)

## End(Not run)

####User-defined target distribution for a calibration problem 
## Not run: 
#Known function which we want to calibrate
CalibFn <- function(y, logit = FALSE) {
  if (logit == TRUE) {
    y <- exp(LogPq(y)$logp)
  }
  out <- 1e3*y[1]^3*y[2]^3/sqrt(20 + y[3])
  return(out)
}

#Generate data 
z <- rnorm(n = 1000, mean = CalibFn(c(1/3, 1/3, 1/3), 2))

#User defined target distribution
Target <- function(ycand, ycurrent, a, dat, pars = NULL) {
  out <- sum(dnorm(dat, CalibFn(ycand, logit = TRUE), 2, log = TRUE)) - 
    sum(dnorm(dat, CalibFn(ycurrent, logit = TRUE), 2, log = TRUE)) + 
    sum((a - 1)*(LogPq(ycand)$logp - LogPq(ycurrent)$logp))
  return(out)
} 

#Run sampler
inputDist <- RunMh(center = c(1/3, 1/3, 1/3), B = 3e4, concentration = 3, 
                   h = c(0.2, 0.2, 0.2), type = 'user', dat = z)

## End(Not run)

[Package SALTSampler version 1.1.0 Index]