# Metro_Hastings: Markov Chain Monte Carlo for Bayesian Inference using... In MHadaptive: General Markov Chain Monte Carlo for Bayesian Inference using adaptive Metropolis-Hastings sampling

## 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

 ```1 2 3``` ```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

## Note

While `Metro_Hastings` has an adaptive proposal structure built in, if `prop_sigma` differs greatly from the covariance structure of the target distribution, stationarity may not be achieved.

## Author(s)

Corey Chivers <corey.chivers@mail.mcgill.ca>

`mcmc_thin`, `plotMH`,`BCI`
 ``` 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 40 41 42 43 44 45``` ```### 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) ```