MHadaptive (version 1.1-8)

Metro_Hastings: Markov Chain Monte Carlo for Bayesian Inference using adaptive Metropolis-Hastings

Description

The function Metro_Hastings performs a general Metropolis-Hastings sampling of a user defined function which returns the un-normalized value (likelihood times prior) of a Bayesian model. The proposal variance-covariance structure is updated adaptively for efficient mixing when the structure of the target distribution is unknown.

Usage

Metro_Hastings(li_func, pars, prop_sigma = NULL,  par_names = NULL, iterations = 50000, burn_in = 1000,  adapt_par = c(100, 20, 0.5, 0.75), quiet = FALSE,...)

Arguments

li_func

user defined function (target distribution) which describes a Bayesian model to be estimated. The function should return the un-normalized log-density function (ie. $log[L(theta | x)P(theta)]$). The first argument to this function should be a vector of parameter values at which to evaluate the function.

pars

vector of initial parameter values defining the starting position of the Markov Chain.

prop_sigma

covariance matrix giving the covariance of the proposal distribution. This matrix need not be positive definite. If the covariance structure of the target distribution is known (approximately), it can be given here. If not given, the diagonal will be estimated via the Fisher information matrix.

par_names

character vector providing the names of each parameter in the model.

iterations

integer: number of iterations to run the chain for. Default 50000.

burn_in

integer: discard the first burn_in values. Default 100.

adapt_par

vector of tuning parameters for the proposal covariance adaptation. Default is c(100, 20, 0.5, 0.75). The first element determines after which iteration to begin adaptation. The second gives the frequency with which updating occurs. The third gives the proportion of the previous states to include when updating (by default 1/2). Finally, the fourth element indicates when to stop adapting (default after 75% of the iterations).

quiet

logical: set to TRUE to suppress printing of chain status.

...

additional arguments to be passed to li_func.

Value

trace

matrix containing the Markov Chain

prop_sigma

adapted covariance matrix of the proposal distribution

par_names

character vector of the parameter names

DIC

Deviance Information Criteria

acceptance_rate

proportion of times proposed jumps were accepted

See Also

mcmc_thin, plotMH,BCI

Examples

### A LINEAR REGRESSION EXAMPLE ####
## Define a Bayesian linear regression model
li_reg<-function(pars,data)
{
    a<-pars[1]      #intercept
    b<-pars[2]      #slope
    sd_e<-pars[3]   #error (residuals)
    if(sd_e<=0){return(NaN)}
    pred <- a + b * data[,1]
    log_likelihood<-sum( dnorm(data[,2],pred,sd_e, log=TRUE) )
    prior<- prior_reg(pars)
    return(log_likelihood + prior)
}

## Define the Prior distributions
prior_reg<-function(pars)
{
    a<-pars[1]          #intercept
    b<-pars[2]          #slope  
    epsilon<-pars[3]    #error

    prior_a<-dnorm(a,0,100,log=TRUE)     ## non-informative (flat) priors on all 
    prior_b<-dnorm(b,0,100,log=TRUE)     ## parameters.  
    prior_epsilon<-dgamma(epsilon,1,1/100,log=TRUE)      

    return(prior_a + prior_b + prior_epsilon)
}

# simulate data
x<-runif(30,5,15)
y<-x+rnorm(30,0,5)
d<-cbind(x,y)


mcmc_r<-Metro_Hastings(li_func=li_reg,pars=c(0,1,1),
    par_names=c('a','b','epsilon'),data=d)

##  For best results, run again with the previously 
##  adapted variance-covariance matrix.

mcmc_r<-Metro_Hastings(li_func=li_reg,pars=c(0,1,1),
    prop_sigma=mcmc_r$prop_sigma,par_names=c('a','b','epsilon'),data=d)

mcmc_r<-mcmc_thin(mcmc_r)
plotMH(mcmc_r)