Here’s some code to animate a rejection sampler

x <- runif(1000)
u <- runif(1000)
hx <- hist(x, plot = FALSE)
M <- 3
post <- x
post[u >= dbeta(x, 2.7, 5.3)/(M*dunif(x))] <- NA
hp <- hist(post, plot = FALSE)
require(rpanel)
redraw <- function(panel){
rp.tkrreplot(panel, tkrp)
panel
}
draw <- function(panel){
r <- as.numeric(panel$Replicates)
par(oma =c(0,0,0,0))
#par(mfrow = c(2,2))
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
curve(dbeta(x, 2.7, 5.3), col = "red", ylab = "f(x) or g(x)", main = "Density Ratio", bty = "n", lwd = 2)
curve(dunif(x, 0, 1), add = TRUE, col = "green", lwd = 2, lty = 1)## was lty 3
arrows(x[r], 0, x[r], dbeta(x[r], 2.7, 5.3), col = "red", lwd = 4, length = 0.01)
arrows(x[r], 0, x[r], dunif(x[r]), col = "green", lwd = 2, lty = 1, length = 0.01) ## was lty 3
##t1 <- bquote(Ratio == .(round(dbeta(x[r], 2.7, 5.3)/(M*dunif(x[r])),2)))
##t1 <- bquote( paste(alpha== .(panel$A), ' ', beta==.(panel$B)))
t1 <- bquote(paste(U == .(round(u[r],2)), '>', Ratio == .(round(dbeta(x[r], 2.7, 5.3)/(M*dunif(x[r])),2))) )
plotcol <- "red"
if (u[r] < dbeta(x[r], 2.7, 5.3)/(M*dunif(x[r]))){
t1 <- bquote(paste(U == .(round(u[r],2)), '<', Ratio == .(round(dbeta(x[r], 2.7, 5.3)/(M*dunif(x[r])),2))) )
plotcol <- "green"
}
text(0.7, 1.9, t1, col = plotcol)
##t2 <- bquote(U == .(u[r]))
##text(0.8, 1.1, t2)
##decision <- "Reject"
##if (u[r] < dbeta(x[r], 2.7, 5.3)/(M*dunif(x[r]))) decision <- "Accept"
##plotcol <- "red"
##if (u[r] < dbeta(x[r], 2.7, 5.3)/(M*dunif(x[r]))) plotcol <- "green"
##text(0.8, 0.7, decision)
points(x[r], 0, pch = 16, col = plotcol)
hist(x[1:r], breaks = hx$breaks, ylim = c(0, max(hx$counts)), main = "Proposal")
points(x[r], 0, pch = 16, col = "green")
hist(post[1:r], xlim = c(0, 1), breaks = hx$breaks, ylim = c(0, max(hp$counts)), main = "Accepted sample")
points(x[r], 0, pch = 16, col = plotcol)
#acplot(1,1)
panel
}
rejectplot <- rp.control(title = "Rejection Sampler", Replicates = 1)
rp.tkrplot(rejectplot, tkrp, draw)
rp.slider(rejectplot, Replicates, action = redraw, from = 1, to = 1000, resolution = 1)

### Like this:

Like Loading...

*Related*