Rizzo, Problem 9.3

Using the unifrom distribution as proposed distribution.

set.seed(8675309)

m <- 10000
x <- numeric(m)
x[1] <- rnorm(1)
k <- 0
u <- runif(m)

for (i in 2:m) {
    xt <- x[i-1]
    y <- rnorm(1, mean = xt)
    num <- dcauchy(y)*dnorm(xt, mean = y)
    den <- dcauchy(xt)*dnorm(y, mean = xt)
    if(u[i] <= num/den) {
        x[i] <- y
    }
    else {
        x[i] <- xt
        k <- k+1
    }
}

Plotting sample after discarding first 1000 of the chain

after_burn <- x[1000:m]
plot(c(1000:m),after_burn, type = "l")

Sampler looks a little unstable at points.

Compare histogram Matro-Hastings Sampler of Cauchy vs standard Cauchy:

s <- seq(0,1,0.01)
t <- qcauchy(s)
t <- t[(t> -Inf) & (t< Inf)]
hist(after_burn, freq = FALSE, main = "", xlab = "x")
lines(t, dcauchy(t), lty = 2, col = "red")

Comparing the deciles:

s <- seq(0,1,.1)
quantile(after_burn, probs = s, type = 5)
##          0%         10%         20%         30%         40%         50% 
## -28.2616799  -3.8625263  -1.5989351  -0.8539173  -0.4628637  -0.1105471 
##         60%         70%         80%         90%        100% 
##   0.2064428   0.5762249   1.0758023   2.1696037  23.4678082
quantile(qcauchy(s), s, type = 5)
##         0%        10%        20%        30%        40%        50% 
##       -Inf       -Inf -1.8867724 -0.8565104 -0.3650820  0.0000000 
##        60%        70%        80%        90%       100% 
##  0.3650820  0.8565104  1.8867724        Inf        Inf

The Matro-Hastings method seems to produce amore of a cluster around the central tendancy, i.e., the tails are shorter than the normal Cauchy