Motivation for Gibbs Sampler

> f.cauchy <- function(x) {
+     p <- (1/pi) * (1/(1 + x^2))
+     return(p)
+ }
> f.normal <- function(x) {
+     sig <- 1
+     mu <- 0
+     d <- sqrt(2 * pi) * sig
+     n <- exp((-1 * (x - mu)^2)/(2 * sig^2))
+     p <- n/d
+     return(p)
+ }
> f.ratio <- function(x) {
+     return(-1 * f.normal(x)/f.cauchy(x))
+ }
> M <- -1 * optim(8, f.ratio)$value
> N <- 10000
> U <- runif(N)
> V <- rcauchy(N, 0)
> normal.var <- ifelse(U < (1/M) * (-1 * f.ratio(V)), V, NA)
> normal.var <- normal.var[!is.na(normal.var)]
> normal.theor <- rnorm(length(normal.var), 0, 1)
> qqplot(normal.var, normal.theor, main = "Check sample vs simulated ")

> qqnorm(normal.var)
> qqline(normal.var, col = "blue", lwd = 2, main = "Check sample vs simulated ")

> f.cauchy <- function(x) {
+     p <- (1/pi) * (1/(1 + x^2))
+     return(p)
+ }
> f.normal <- function(x) {
+     sig <- 1
+     mu <- 0
+     d <- sqrt(2 * pi) * sig
+     n <- exp((-1 * (x - mu)^2)/(2 * sig^2))
+     p <- n/d
+     return(p)
+ }
> f.ratio <- function(x) {
+     return(-1 * f.cauchy(x)/f.normal(x))
+ }
> M <- -1 * optim(10, f.ratio)$value
> N <- 10000
> U <- runif(N)
> V <- rnorm(N, 0, 1)
> cauchy.var <- ifelse(U < (1/M) * (f.ratio(V)), V, NA)
> cauchy.var <- cauchy.var[!is.na(cauchy.var)]