a = 3
b = 10
par(mfrow=c(2,2))
for(m in c(10,20,100,1000))
{
c = 1/m
xdom = c(0:(m-1))/m +c/2
x_start = xdom - c/2
x_end = xdom + c/2
prob_x = dbeta(xdom ,a,b)
prob_x = prob_x/sum(prob_x)
altezza = prob_x/c
xseq = seq(0,1, by = 0.001)
plot(xseq, dbeta(xseq,a,b), type="l", main=paste("m=", m), lwd=1)
rect(x_start, rep(0, m),x_end,altezza , border=2, col=5, lwd= c*10)
points(xdom, rep(0, m), pch=20, cex= c*5, col=2)
lines(xseq, dbeta(xseq,a,b), col=1, lwd=1)
}
par(mfrow=c(1,1))