SIR.for.epsilon.1 <- function(a1, b1, a2, b2, x1, n1, x2, n2, size=10000) { theta1 <- runif(size, min=0,max=1) epsilon <- runif(size, min=theta1-1, max=theta1) weight <- theta1^(a1+x1-1)*(1-theta1)^(b1+n1-x1-1)* ( theta1 - epsilon)^(a2+x2-1)*(1- (theta1 - epsilon))^(b2+n2-x2-1) epsilon.post <- sample(epsilon, size, replace=T, prob=weight) layout(matrix(c(1,2), byrow=T, nrow=2)) hist(epsilon.post) hist(weight) return(quantile(epsilon.post, prob=c(0.025, 0.25, 0.5, 0.75, 0.975))) }