Rejection Sampler

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)

 

Advertisements
This entry was posted in Uncategorized. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s