# 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

1 |

### Arguments

`center` |
Vector of numeric values summing to 1 that define the center of the distributional parameters of the posterior. For type |

`B` |
Number of iterations to run the chain |

`concentration` |
This argument specifies the concentration parameter where |

`h` |
Vector of step sizes. Length of vector must match length of |

`type` |
Specifies the target distribution. Select type |

`dat` |
A matrix or vector passing data to the sampler. For type |

`pars` |
A list of additional parameters that can be passed to the user-specified target function for type |

### 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, *θ* = `ycand`

, to the log-likelihood of the posterior for the current value, *θ* = `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 of numeric values summing to |

`B` |
Number of iterations to run the chain |

`concentration` |
For type |

`h` |
Vector of step sizes. Length of vector must match length of |

`type` |
Specifies the target distribution. Select type |

`dat` |
A matrix or vector passing data to the sampler. For type |

`a` |
Dirichlet distribution parameters, |

### Examples

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 | ```
###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)
``` |