Week10

---

Rizzo 9.3

Use the Metropolis-Hastings sampler to generate random variables from a standard Cauchy distribution. Discard the first 1000 of the chain, and compare 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( \theta, \eta ) \] distribution has density function

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

---

Generate random variables from a standard Cauchy distribution

num_iter <- 5000
my_vector <- numeric(num_iter)
my_vector[1] <- rnorm(1)
counter <- 0
#generates random deviates
u <- runif(num_iter)

for (i in 2:num_iter) 
{
    xt <- my_vector[i-1]
    y <- rnorm(1, mean = xt)
    #Density Cauchy distribution
    numerator <- dcauchy(y)*dnorm(xt, mean = y)
    denominator <- dcauchy(xt)*dnorm(y, mean = xt)
    
    if(u[i] <= numerator/denominator) 
    {
        my_vector[i] <- y
    }
    else 
    {
        my_vector[i] <- xt
        counter <- counter+1
    }
}

plot(my_vector,type = "l")

after1kDiscard <- my_vector[1000:num_iter]

Sampler(after discarding1k ) vs standard Cauchy

#Sequence Generation
sequence <- seq(0,1,0.01)
standardCauchy <- qcauchy(sequence)
standardCauchy <- standardCauchy[(standardCauchy> -Inf) & (standardCauchy< Inf)]
hist(after1kDiscard, freq = FALSE,main = "")
lines(standardCauchy, dcauchy(standardCauchy), lty = 2)