N = 9; x = 0:N; s1=2; s2=3
dy = dbetabin.ab(x, size=N, shape1=s1, shape2=s2)
barplot(rbind(dy, dbinom(x, size=N, prob=s1/(s1+s2))),
beside = TRUE, col = c("blue","green"),
main=paste("Beta-binomial (size=",N,", shape1=",s1,
", shape2=",s2,") (blue) vs\n",
" Binomial(size=", N, ", prob=", s1/(s1+s2), ") (green)", sep=""),
names.arg = as.character(x))
sum(dy*x) # Check expected values are equal
sum(dbinom(x, size=N, prob=s1/(s1+s2))*x)
cumsum(dy) - pbetabin.ab(x, N, shape1=s1, shape2=s2)
y = rbetabin.ab(n=10000, size=N, shape1=s1, shape2=s2)
ty = table(y)
barplot(rbind(dy, ty/sum(ty)),
beside = TRUE, col = c("blue","red"),
main=paste("Beta-binomial (size=",N,", shape1=",s1,
", shape2=",s2,") (blue) vs\n",
" Random generated beta-binomial(size=", N, ", prob=", s1/(s1+s2),
") (red)", sep=""),
names.arg = as.character(x))
Run the code above in your browser using DataLab