set.seed(1); rbetabinom(10, 100, prob = 0.5)
set.seed(1); rbinom(10, 100, prob = 0.5) # The same since rho = 0
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 = 10000, 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))
set.seed(208); N <- 1000000; size = 20;
pstr0 <- 0.2; pstrsize <- 0.2
k <- rozibetabinom.ab(N, size, s1, s2, pstr0, pstrsize)
hist(k, probability = TRUE, border = "blue",
main = "Blue = inflated; orange = ordinary beta-binomial",
breaks = -0.5 : (size + 0.5))
sum(k == 0) / N # Proportion of 0
sum(k == size) / N # Proportion of size
lines(0 : size,
dbetabinom.ab(0 : size, size, s1, s2), col = "orange")
lines(0 : size, col = "blue",
dozibetabinom.ab(0 : size, size, s1, s2, pstr0, pstrsize))
Run the code above in your browser using DataLab