Rizzo, Maria L. 2008. Statistical Computing with R. Chapter 9: Problem 3.

Cauchy(\(\theta, \eta\)) distribution has PDF:

\(f(x) = \frac{1}{\theta \pi (1 + [(x - \eta)/ \theta]^2)}, -\infty < x < \infty, \theta > 0\)

Standard Cauchy has the Cauchy(\(\theta = 1, \eta = 0\)) density

Metropolis-Hastings sampler for standard Cauchy distribution:

Using proposal distribution N(\(\mu = X_t, \sigma^2 = 1\))

library(ggplot2)

f <- function(x, theta, eta) {
  stopifnot(theta > 0)
  d <- 1 / (theta * pi * (1 + ((x - eta) / theta)^2))
  return(d)
}

m <- 10000
theta = 1
eta = 0
x <- numeric(m)
x[1] <- rnorm(1, mean = 0, sd = 1)
k <- 0 
u <- runif(m)

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

x <- x[-(1:1000)]  # remove first 1000 from Markov chain

ggplot(data = data.frame(x = x), aes(x = seq_along(x), y = x)) +
  geom_line() +
  ylab("index")

ggplot(data = data.frame(x = x), aes(x)) + 
  geom_histogram(aes(y = ..density..), bins = 100) +
  stat_function(fun = dcauchy, args = list(scale = theta, location = eta), color = "red")

ob <- quantile(x, seq(0, 1, 0.1))
ex <- qcauchy(seq(0, 1, 0.1), scale = theta, location = eta)
df <- data.frame(decile = names(ob), observed = ob, expected = ex, row.names = NULL)
knitr::kable(df)
decile observed expected
0% -14.2974037 -Inf
10% -2.1395533 -3.0776835
20% -1.1210559 -1.3763819
30% -0.5817413 -0.7265425
40% -0.2340947 -0.3249197
50% 0.0648931 0.0000000
60% 0.3642279 0.3249197
70% 0.7341899 0.7265425
80% 1.3615419 1.3763819
90% 2.7835694 3.0776835
100% 13.2769199 Inf