rm(list=ls()); set.seed(20210127);
# Data
n1 <- 10; y1 <- 5 # Data from experiment 1 n2 <- 10; y2 <- 8 # Data from experiment 2 S <- 1000 # Number of MC samples theta1 <- rbeta(S, y1+1, n1-y1+1) theta2 <- rbeta(S, y2+1, n2-y2+1) # Compare histograms for posteriors of theta1 and theta2 op <- par(mfrow=c(1,2)) hist(theta1, freq=F, right=F, xlab=expression(theta[1]), xlim=0:1, breaks=25, main="MC approximate posterior") hist(theta2,, freq=F, right=F, xlab=expression(theta[2]), xlim=0:1, breaks=25, main="MC approximate posterior") par(op) # Scatterplot of Monte Carlo sample # Can use to approximate posterior prob that theta2>theta1
(y1+1)/(n1+2) # True post mean
mean(theta1) # MC estimate
plot(theta1, theta2, xlim=c(0,1), ylim=c(0,1),
xlab=expression(theta[1]), ylab=expression(theta[2]))
abline(0,1)
mean(theta2>theta1)