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)
}