# NOT RUN {
set.seed(1); rbetabinom(10, 100, prob = 0.5)
set.seed(1); rbinom(10, 100, prob = 0.5) # The same since rho = 0
# }
# NOT RUN {
N <- 9; xx <- 0:N; s1 <- 2; s2 <- 3
dy <- dbetabinom.ab(xx, size = N, shape1 = s1, shape2 = s2)
barplot(rbind(dy, dbinom(xx, size = N, prob = s1 / (s1+s2))),
beside = TRUE, col = c("blue","green"), las = 1,
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(xx), cex.main = 0.8)
sum(dy * xx) # Check expected values are equal
sum(dbinom(xx, size = N, prob = s1 / (s1+s2)) * xx)
cumsum(dy) - pbetabinom.ab(xx, N, shape1 = s1, shape2 = s2) # Should be all 0
y <- rbetabinom.ab(n = 1e4, size = N, shape1 = s1, shape2 = s2)
ty <- table(y)
barplot(rbind(dy, ty / sum(ty)),
beside = TRUE, col = c("blue", "orange"), las = 1,
main = paste("Beta-binomial (size=", N, ", shape1=", s1,
", shape2=", s2, ") (blue) vs\n",
" Random generated beta-binomial(size=", N, ", prob=", s1/(s1+s2),
") (orange)", sep = ""), cex.main = 0.8,
names.arg = as.character(xx))
N <- 1e5; size <- 20; pstr0 <- 0.2; pstrsize <- 0.2
kk <- rzoibetabinom.ab(N, size, s1, s2, pstr0, pstrsize)
hist(kk, probability = TRUE, border = "blue", ylim = c(0, 0.25),
main = "Blue/green = inflated; orange = ordinary beta-binomial",
breaks = -0.5 : (size + 0.5))
sum(kk == 0) / N # Proportion of 0
sum(kk == size) / N # Proportion of size
lines(0 : size,
dbetabinom.ab(0 : size, size, s1, s2), col = "orange")
lines(0 : size, col = "green", type = "b",
dzoibetabinom.ab(0 : size, size, s1, s2, pstr0, pstrsize))
# }
Run the code above in your browser using DataLab