Discussion 10

Use the Metropolis-Hastings sampler to generate random variables from a standard Cauchy distribution. Discard the first 1000 of the chain, and com- pare the deciles of the generated observations with the deciles of the standard Cauchy distribution (see qcauchy or qt with df=1). Recall that a Cauchy(θ, η) distribution has density function

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

The standard Cauchy has the Cauchy(θ = 1, η = 0) density. (Note that the standard Cauchy density is equal to the Student t density with one degree of freedom.)

set.seed(2017)
f <- function(x, sigma, eta) {
if (any(x < 0)) return (0)
stopifnot(sigma > 0)
return(1/(sigma * pi * (1 + ((x - eta)/sigma)^2) ))
}


m <- 3000
eta <- 0
sigma <- 1
x <- numeric(m)
x[1] <- rchisq(1, df=1)
k <- 0
u <- runif(m)

for (i in 2:m) {
  xt <- x[i-1]
  y <- rchisq(1, df=xt)
  num <- f(y, sigma, eta) * dchisq(xt, df=y)
  den <- f(xt, sigma, eta) * dchisq(y, df=xt)
  if (u[i] <= num/den) x[i] <- y else {
    x[i] <- xt
    k <- k+1
  }
}
xL <- x[1001:length(x)]

summary(x) #before removing first 1k

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.03032  0.45364  0.96936  2.27926  2.09353 57.08266

summary(xL)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.03032  0.52721  1.08846  2.71790  2.41876 57.08266

summary(qcauchy(runif(1000)))

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## -2244.6049    -1.0234    -0.0875    -3.5680     0.9157   181.2768

#Some similarities?

#Plots:
plot(density(xL))

plot(density(dcauchy(runif(1000))))