Problem 9.3 from Statistical Computing with R: Rizzo, Maria
Metropolis-Hastings sampler for standard Cauchy Distribution
set.seed(123)
m <- 5000
x <- numeric(m)
x[1] <- rnorm(1)
rejected <- 0
#generates random deviates
u <- runif(m)
for (i in 2:m) {
xt <- x[i-1]
y <- rnorm(1, mean = xt)
#Density Cauchy distribution
num <- dcauchy(y) * dnorm(xt, mean = y)
den <- dcauchy(xt) * dnorm(y, mean = xt)
if(u[i] <= num/den) {
x[i] <- y
}
else
{
x[i] <- xt
rejected <- rejected + 1
}
}
Rejections
- Nubmer rejected = 1219
- Total Samples = 5000
- % Reject = 24.38
Plots
# discard the first 1000
y1 <- x[1001:m]
plot(y1 ,type = "l", main="", ylab="x")
sequence <- seq(0,1,0.01)
cauchy <- qcauchy(sequence)
hist(y1, main="", freq = FALSE)
lines(cauchy, dcauchy(cauchy), lty = 2)
Deciles - Generated Observations vs. Standard Cauchy Distribution
observed <- quantile(y1, seq(0, 1, 0.1))
standard <- qcauchy(seq(0, 1, 0.1))
knitr::kable(cbind(observed, standard))
| 0% |
-10.8073155 |
-Inf |
| 10% |
-2.2483208 |
-3.0776835 |
| 20% |
-1.2073383 |
-1.3763819 |
| 30% |
-0.6690567 |
-0.7265425 |
| 40% |
-0.3116071 |
-0.3249197 |
| 50% |
-0.0311638 |
0.0000000 |
| 60% |
0.2830026 |
0.3249197 |
| 70% |
0.6694872 |
0.7265425 |
| 80% |
1.2600980 |
1.3763819 |
| 90% |
2.8869523 |
3.0776835 |
| 100% |
12.5638759 |
Inf |
