MH = function(start,hfun,Jfun,canfun,niter) { ## Function for a Metropolis-Hastings algorithm ## The arguments to the function: ## start - is a vector of starting values for the parameters ## hfun - is a function for computing the log(h(X)) function in M-H, ## usually this is a likelihood or posterior. ## Note that this should be the log of the likelihood/posterior ## The arguments to this function should be the parameters ## There should be a single return value ## Jfun - is the jumping distribution ## Note that this should be the log of the Jumping distribution ## The arguments to this function should be the parameters, even though ## they may not actually be used ## There should be a single return value, return a constant value ## if using a Metropolis algorithm ## canfun - a function to generate candidate values, it should return a vector ## of candidate values corresponding to start ## there should be no arguments to this function, only return values ## niter - number of iterations for the algorithm ## ## The function returns an n by p matrix, where n = niter and p = number of ## parameters # Set up a matrix of parameter values at each iteration parsi = matrix(0,nrow=niter,ncol=length(start)) cpars = start # Loop through iterations of M-H for (i in 1:niter) { # Get candidate values candvals = canfun() # Compute the log of the acceptance ratio r logr = hfun(candvals) - Jfun(candvals) - hfun(cpars) + Jfun(cpars) # Generate probability and check for acceptance logu = log(runif(1)) if (logu < logr) { cpars = candvals } parsi[i,] = cpars } return(parsi) }