Cuny 604 - Discussion 10

Problem 9.3 page 277

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(??, ??) distribution has density function of:

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

The choice of proposal distribution is very flexible, but the chain generated by this choice must satisfy certain regularity conditions. The proposal distribution must be chosen so that the generated chain will converge to a stationary distribution - the target distribution f. Required conditions for the generated chain are irreducibility, positive recurrence, and aperiodicity

1. Choose a proposal distribution g(·|Xt) (subject to regularity conditions stated above), in this case, normal distribution
   
2. Generate X0 from a distribution g.
   
3. Repeat (until the chain has converged to a stationary distribution according to some criterion):
   

• Generate Y from g(.Xt)
  
• Generate U from [0,1]
  
• \(if\quad U\le \frac { f\left( Y \right) g(X_{ t }|Y) }{ f({ X }_{ t })g(Y|{ X }_{ t }) } ,\quad accept\quad Y\quad and\quad set\quad { X }_{ t+1 }\quad to\quad Y,\quad otherwise\quad { X }_{ t+1 }={ X }_{ t }\)
  

4. increment t

Nsim <- 10^4  
X <- c(rt(1,1)) # initialize the chain from the stationary 

for (t in 2:Nsim){
  Y <-  rnorm(1) # candidate normal 
  U <-runif(1)
  rho <- dt(Y,1)*dnorm(X[t-1])/(dt(X[t-1],1)*dnorm(Y)) 
  X[t]<- X[t-1] + (Y-X[t-1])*(U<=rho)  #U<=rho either 0 or 1
}

# Discard first 1000 value of chain 
X_keep <- X[10001:length(X)]

Estimated_values <- quantile(X, seq(0, 1, 0.1))
Actual_values <- qcauchy(seq(0, 1, 0.1))

results <- cbind.data.frame(Estimated_values, Actual_values)

results

##      Estimated_values Actual_values
## 0%        -3.57949284          -Inf
## 10%       -1.62055266    -3.0776835
## 20%       -0.94275139    -1.3763819
## 30%       -0.52764746    -0.7265425
## 40%       -0.22784786    -0.3249197
## 50%        0.03510771     0.0000000
## 60%        0.30067721     0.3249197
## 70%        0.61990386     0.7265425
## 80%        1.05790353     1.3763819
## 90%        1.91088093     3.0776835
## 100%       3.52360467           Inf

References:
Rizzo; Maria L.. Statistical Computing with R (Chapman & Hall/CRC The R Series)
www.springer.com/cda/content/document/cda…/9781441915757-c1.pdf?SGWID