mh.fit: Fitter function for Metropolis-Hastings (MH)

In sthomas522/hmclearn: Fit Statistical Models Using Hamiltonian Monte Carlo

Description

This is the basic computing function for MH and should not be called directly except by experienced users.

Usage

mh.fit(
  N,
  theta.init,
  qPROP,
  qFUN,
  logPOSTERIOR,
  nu = 0.001,
  varnames = NULL,
  param = list(),
  ...
)

Arguments

N	Number of MCMC samples
theta.init	Vector of initial values for the parameters
qPROP	Function to generate proposal
qFUN	Probability for proposal function. First argument is where to evaluate, and second argument is the conditional parameter
logPOSTERIOR	Function to calculate and return the log posterior given a vector of values of theta
nu	Single value or vector parameter passed to qPROP or qFUN for the proposal density
varnames	Optional vector of theta parameter names
param	List of additional parameters for logPOSTERIOR
...	Additional parameters for logPOSTERIOR

Value

List for mh

Elements in mh list

N

Number of MCMC samples

theta

Nested list of length N of the sampled values of theta for each chain

thetaCombined

List of dataframes containing sampled values, one for each chain

r

NULL for Metropolis-Hastings

theta.all

Nested list of all parameter values of theta sampled prior to accept/reject step for each

r.all

NULL for Metropolis-Hastings

accept

Number of accepted proposals. The ratio accept / N is the acceptance rate

accept_v

Vector of length N indicating which samples were accepted

M

NULL for Metropolis-Hastings

algorithm

MH for Metropolis-Hastings

Examples

# Logistic regression example
X <- cbind(1, seq(-100, 100, by=0.25))
betavals <- c(-0.9, 0.2)
lodds <- X %*% betavals
prob1 <- as.numeric(1 / (1 + exp(-lodds)))

set.seed(9874)
y <- sapply(prob1, function(xx) {
  sample(c(0, 1), 1, prob=c(1-xx, xx))
})

f1 <- mh.fit(N = 2000,
         theta.init = rep(0, 2),
         nu = c(0.03, 0.001),
         qPROP = qprop,
         qFUN = qfun,
         logPOSTERIOR = logistic_posterior,
         varnames = paste0("beta", 0:1),
         y=y, X=X)

f1$accept / f1$N

sthomas522/hmclearn documentation built on Oct. 16, 2020, 8:34 a.m.