Example of a Cauchy Gibbs Sampler

GibbsSamplierExample.R

library(MCMCpack)
rm(list=ls())
set.seed(83247279)
reps=50000

# draw data from the cauchy distribution
x.dat <- rcauchy(100, location=2, scale=0.8)


######### now recover the location and the scale


# likelihood: the cauchy distribution

cauchy <- function(x, x0, gamma){
  
  if (gamma <= 0) {
    return(0)
  }
  
  out <- 1/(pi*gamma*(1+((x-x0)/gamma)^2))
  
  return(prod(out))
  
}

# prior over the location

prior <- function(x0, gamma){
  
  out<-c()
  out[1] <- dnorm(x0, mean=0, sd=3)
  
  if (gamma > 0) {
    out[2] <- dinvgamma(gamma, shape=1e-6, scale=1)
  } else {
    out[2] <- 0
  }
  
  return(prod(out))
  
}


# the posterior function

cauchy.post <- function(x0, gamma, x){
  out <- prior(x0=x0, gamma=gamma)*cauchy(x=x, x0=x0, gamma=gamma)
  return(out)
}


chain.1 <- c(0)
chain.2 <- c(1)

for(i in 1:reps){
  
  proposal <- chain.1[i] + runif(1, min=-2, max=2)
  num <- cauchy.post(x0=proposal, gamma=chain.2[i], x=x.dat)
  den <- cauchy.post(x0=chain.1[i], gamma=chain.2[i], x=x.dat)
  accept <- runif(1) < num/den
  chain.1[i+1] <- ifelse(accept==T, proposal, chain.1[i])
  
  proposal <- chain.2[i] + runif(1, min=-2, max=2)
  num <- cauchy.post(x0=chain.1[i+1], gamma=proposal, x=x.dat)
  den <- cauchy.post(x0=chain.1[i+1], gamma=chain.2[i], x=x.dat)
  accept <- runif(1) < num/den
  chain.2[i+1] <- ifelse(accept==T, proposal, chain.2[i])
  
}

plot(density(chain.1[1000:50000]), main="x0")
plot(density(chain.2[1000:50000]), main="gamma")

quantile(chain.1[1000:50000], probs=(c(0.025, 0.975)))
quantile(chain.2[1000:50000], probs=(c(0.025, 0.975)))

https://www.youtube.com/watch?v=j4nEAqUUnVw