Fit a generic model using Metropolis-Hastings (MH)

This function runs the MH algorithm on a generic model provided the logPOSTERIOR function. All parameters specified within the list param are passed to these the posterior function.

mh(
  N,
  theta.init,
  qPROP,
  qFUN,
  logPOSTERIOR,
  nu = 0.001,
  varnames = NULL,
  param = list(),
  chains = 1,
  parallel = FALSE,
  ...
)

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 and glogPOSTERIOR
chains	Number of MCMC chains to run
parallel	Logical to set whether multiple MCMC chains should be run in parallel
...	Additional parameters for logPOSTERIOR

Value

Object of class hmclearn

Elements for hmclearn objects

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

varnames

Optional vector of parameter names

chains

Number of MCMC chains

Available logPOSTERIOR functions

linear_posterior

Linear regression: log posterior

logistic_posterior

Logistic regression: log posterior

poisson_posterior

Poisson (count) regression: log posterior

lmm_posterior

Linear mixed effects model: log posterior

glmm_bin_posterior

Logistic mixed effects model: log posterior

glmm_poisson_posterior

Poisson mixed effects model: log posterior

Author

Samuel Thomas samthoma@iu.edu, Wanzhu Tu wtu@iu.edu

Examples

# Linear regression example
set.seed(521)
X <- cbind(1, matrix(rnorm(300), ncol=3))
betavals <- c(0.5, -1, 2, -3)
y <- X%*%betavals + rnorm(100, sd=.2)

f1_mh <- mh(N = 3e3,
         theta.init = c(rep(0, 4), 1),
         nu <- c(rep(0.001, 4), 0.1),
         qPROP = qprop,
         qFUN = qfun,
         logPOSTERIOR = linear_posterior,
         varnames = c(paste0("beta", 0:3), "log_sigma_sq"),
         param=list(y=y, X=X), parallel=FALSE, chains=1)

summary(f1_mh, burnin=1000)
#> Summary of MCMC simulation
#> 
#>                    2.5%         5%        25%        50%        75%        95%
#> beta0         0.4927326  0.5002949  0.5233041  0.5355716  0.5491311  0.5749760
#> beta1        -1.0427050 -1.0366202 -1.0210319 -1.0069922 -0.9976423 -0.9784314
#> beta2         1.9825087  1.9846028  2.0037900  2.0183856  2.0312183  2.0589282
#> beta3        -3.0259491 -3.0182092 -3.0008568 -2.9824039 -2.9648561 -2.9412003
#> log_sigma_sq -3.2972533 -3.2437557 -3.1452192 -3.0795228 -2.9635472 -2.7779556
#>                   97.5%
#> beta0         0.5763848
#> beta1        -0.9715624
#> beta2         2.0653425
#> beta3        -2.9338422
#> log_sigma_sq -2.7132118