## A likelihood that arises from observations of 3 phenotypes with ## probabilities c(2+theta, 1-theta, theta)/4, hence a posterior density ## for a uniform prior. f <- function(theta) (2+theta)^14*(1-theta)^3*theta^5 slice_sampler <- function(f, M, L = 0, R = 1, plotit = TRUE, hist = TRUE, delay = 0.1) { x <- seq(L, R, length = 1000) area <- sum(f(x)) * (R-L)/1000 if (plotit) plot(x, f(x)/area, xlim=c(L,R), type = "l", xlab="", ylab="density") res <- numeric(M) X <- runif(1, L, R); Y <- 0 points(X, Y) for (i in 1:M) { Yprev <- Y Y <- runif(1, 0, f(X)/area) if(plotit) { lines(c(X,X), c(Yprev, Y)) points(X, Y) } Xprev <- X repeat { # use unimodality X <- runif(1, L, R) if (f(X) > Y * area) break } res[i] <- X if(plotit) { lines(c(Xprev, X), c(Y, Y)) points(X, Y) Sys.sleep(delay) } } res } samp <- slice_sampler(f, 10, delay=1) samp <- slice_sampler(f, 10000, plotit=FALSE) L <- 0; R <- 1 x <- seq(L, R, length = 1000) area <- sum(f(x)) * (R-L)/1000 hist(samp, prob = TRUE, breaks = seq(0, 1, by=0.01)) lines(x, f(x)/area)