MH algorithm to generate beta random variables

[GaryBAYLOR]

We compare the density of random samples by using MH algorithm and that by using base R function rbeta. The code is as follows.

mh.beta <- function(a = 2.7, b = 6.3, N = 1e4) {
    X <- numeric(N)
    X[1] <- runif(1)

    for(i in 2:N) {
        y <- runif(1)
        rho <- dbeta(y, a, b) / dbeta(X[i-1], a, b)
        X[i] <- X[i-1] + (y - X[i-1]) * (runif(1) < rho)
    }

    plot(density(X), type = "l", main = "MH simulation & Direct simulation Comparision", xlab = "")
    lines(density(rbeta(N, a, b)), type = "l", col = "red")
    legend(0.65, 2.65, legend = c("MH simulation", "Direct simulation"), lty = c(1, 1), col = c("black", "red"))
}

The plot when N = 1e5 is like the following.

We can see that MH is performing quite good, except near the peak.