Problem 9.3 Metropolis Hastings Sampler

formula

formula

#Sample size
N=10^4

x <- numeric(N) 
#generating the first value from the normal distribution.
x[1] <- rnorm(1) 
#generate U from uniform[0,1]
u <- runif(N) 
#Intilialize the rejects
k <- 0 

for (i in 2:N) { 
  
  y <- rnorm(1, x[i-1]) 
  #Now using the density for cauchy distribution.
  if (u[i] <= dcauchy(y) * dnorm(x[i-1],mean = y)/dcauchy(x[i-1]) * dnorm(y,mean = x[i-1]))   {
    x[i] <- y 
  }
  else { 
    x[i] <- x[i-1] 
    #incrementing the rejects
    k <- k + 1
  }
} 

k#the value of rejects
## [1] 9069
x=x[1001:length(N)]


# Qns: Comparing the denciles of generated Observations with standard cauchy distribution
obseved = quantile(x,seq(0,1,0.1))
expected = qcauchy(seq(0,1,0.1))
knitr::kable(data.frame(obseved,expected
                ))
obseved expected
0% -1.2456778 -Inf
10% -0.6938841 -3.0776835
20% -0.2551065 -1.3763819
30% -0.1319888 -0.7265425
40% 0.0214121 -0.3249197
50% 0.2009202 0.0000000
60% 0.2481482 0.3249197
70% 0.5888677 0.7265425
80% 0.7690502 1.3763819
90% 1.7115989 3.0776835
100% 3.3201444 Inf