R/Metropolis.R

In Scopia/StatComp18053: generating and monitoring convergence of Markov chain and permutation test for equal distribution

#' @title generate Markov chains for specific distribution based on Normal distribution
#' @description This function generates a Markov chain using a random walk Metropolis-Hastings algorithm.The user supplies target distribution, burn-in time, the length of the chain and the variance of proposal distribution. And a Markov chain after discarding burn-in samples is returned, which can be used for monitoring convengence.
#' @usage Metropolis(burn_in=1000,dist=dcauchy,sigma,N=10000,print_acc=F)
#' @param burn_in the total length of discarding
#' @param dist the target distribution
#' @param sigma the variance of the normal proposal distribution
#' @param N the length of Markov chain
#' @param print_acc if print acceptance rate or not
#' @details Metropolis generates a Markov chain using a Metropolis-Hastings algorithm. The aim is to generate random numbers from specific distribution based on Normal proposal distribution
#' @references Statistic computing with R. Maria L. Rizzo
#' @examples
#' \dontrun{
#' y=Metropolis(sigma=2,print_acc=T)
#' plot(density(Metropolis(sigma = 3,print_acc = T)))
#' }
#' @export
Metropolis=function(burn_in=1000,dist=dcauchy,sigma,N=10000,print_acc=F){
  y=1:N
  y[1] <- rnorm(1,sd=sigma)
  d=dist
  acc=0
  for(i in 2:N){
    yt=y[i-1]
    x <- rnorm(1, mean = yt,sd=sigma)
    num <- d(x) * dnorm(yt, mean = x,sd=sigma)
    den <- d(yt) * dnorm(x, mean = yt,sd=sigma)
    if(runif(1)<=min(abs(num/den),1)){y[i]=x;acc=acc+1}
    else y[i]=yt
  }
  if(print_acc)
  {print(acc/N)}
  return(y[(burn_in+1):N])
}