DATA 604 Discussion 10

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

N <- 5000
x <- numeric(N)
x[1] <- rnorm(1)
k <- 0
u <- runif(N)
for (i in 2:N){
  y <- rnorm(1, mean = x[i-1])
  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]
      k <- k + 1
      }
  }
  
x <- x[1001:N]
plot(x,type = "l")

observed <- quantile(x, seq(0, 1, 0.1))
expected <- qcauchy(seq(0, 1, 0.1))
c<-cbind(observed, expected)
c

##        observed   expected
## 0%   -2.3063389       -Inf
## 10%  -0.7200077 -3.0776835
## 20%  -0.4677933 -1.3763819
## 30%  -0.1924404 -0.7265425
## 40%  -0.0352437 -0.3249197
## 50%   0.1730855  0.0000000
## 60%   0.3233978  0.3249197
## 70%   0.5305910  0.7265425
## 80%   0.7675695  1.3763819
## 90%   1.3952413  3.0776835
## 100%  4.1123138        Inf