SIR.for.epsilon.4 <- function(a1, b1, a2, b2, x1, n1, x2, n2, size=10000) { theta1 <- rbeta(size, a1+x1, b1+n1-x1) epsilon <- rbeta(size, 30,70) weight <- theta1^(a1+x1-1)*(1-theta1)^(b1+n1-x1-1)* ( theta1 - epsilon)^(a2+x2-1)* (1- (theta1 - epsilon))^(b2+n2-x2-1)/(dbeta(epsilon, 30, 70) *dbeta(theta1, a1+x1, b1+n1-x1)) epsilon.post<-sample(epsilon, length(epsilon), 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))) }