The Beta-Binomial Distribution

Density, distribution function, and random generation for the beta-binomial distribution and the inflated beta-binomial distribution.

Usage

dbetabinom(x, size, prob, rho = 0, log = FALSE)
pbetabinom(q, size, prob, rho = 0, log.p = FALSE)
rbetabinom(n, size, prob, rho = 0)
dbetabinom.ab(x, size, shape1, shape2, log = FALSE,
              Inf.shape = exp(20), limit.prob = 0.5)
pbetabinom.ab(q, size, shape1, shape2, limit.prob = 0.5,
              log.p = FALSE)
rbetabinom.ab(n, size, shape1, shape2, limit.prob = 0.5,
              .dontuse.prob = NULL)
dzoibetabinom(x, size, prob, rho = 0, pstr0 = 0, pstrsize = 0,
              log = FALSE)
pzoibetabinom(q, size, prob, rho, pstr0 = 0, pstrsize = 0,
              lower.tail = TRUE, log.p = FALSE)
rzoibetabinom(n, size, prob, rho = 0, pstr0 = 0, pstrsize = 0)
dzoibetabinom.ab(x, size, shape1, shape2, pstr0 = 0, pstrsize = 0,
                 log = FALSE)
pzoibetabinom.ab(q, size, shape1, shape2, pstr0 = 0, pstrsize = 0,
              lower.tail = TRUE, log.p = FALSE)
rzoibetabinom.ab(n, size, shape1, shape2, pstr0 = 0, pstrsize = 0)

Arguments

x, q: vector of quantiles.

size

number of trials.

n

number of observations. Same as runif.

prob

the probability of success $\mu$. Must be in the unit closed interval $[0,1]$.

rho

the correlation parameter $\rho$, which should be in the interval $[0, 1)$. The default value of 0 corresponds to the usual binomial distribution with probability prob. Setting rho = 1 would set both shape parameters equal to 0, and the ratio 0/0, which is actually NaN, is interpreted by Beta as 0.5. See the warning below.

shape1, shape2

the two (positive) shape parameters of the standard beta distribution. They are called a and b in beta respectively. Note that shape1 = prob*(1-rho)/rho and shape2 = (1-prob)*(1-rho)/rho is an important relationship between the parameters, so that the shape parameters are infinite by default because rho = 0; hence limit.prob = prob is used to obtain the behaviour of the usual binomial distribution.

log, log.p, lower.tail

Same meaning as runif.

Inf.shape

Numeric. A large value such that, if shape1 or shape2 exceeds this, then special measures are taken, e.g., calling dbinom. Also, if shape1 or shape2 is less than its reciprocal, then special measures are also taken. This feature/approximation is needed to avoid numerical problem with catastrophic cancellation of multiple lbeta calls.

limit.prob

Numerical vector; recycled if necessary. If either shape parameters are Inf then the binomial limit is taken, with shape1 / (shape1 + shape2) as the probability of success. In the case where both are Inf this probability will be a NaN = Inf/Inf, however, the value limit.prob is used instead. Hence the default for dbetabinom.ab() is to assume that both shape parameters are equal as the limit is taken (indeed, Beta uses 0.5). Note that for [dpr]betabinom(), because rho = 0 by default, then limit.prob = prob so that the beta-binomial distribution behaves like the ordinary binomial distribution with respect to arguments size and prob.

.dontuse.prob

An argument that should be ignored and not used.

pstr0

Probability of a structual zero (i.e., ignoring the beta-binomial distribution). The default value of pstr0 corresponds to the response having a beta-binomial distribuion inflated only at size.

pstrsize

Probability of a structual maximum value size. The default value of pstrsize corresponds to the response having a beta-binomial distribution inflated only at 0.

Value

dbetabinom and dbetabinom.ab give the density, pbetabinom and pbetabinom.ab give the distribution function, and rbetabinom and rbetabinom.ab generate random deviates.

dzoibetabinom and dzoibetabinom.ab give the inflated density, pzoibetabinom and pzoibetabinom.ab give the inflated distribution function, and rzoibetabinom and rzoibetabinom.ab generate random inflated deviates.

Author

T. W. Yee and Xiangjie Xue

Details

The beta-binomial distribution is a binomial distribution whose probability of success is not a constant but it is generated from a beta distribution with parameters shape1 and shape2. Note that the mean of this beta distribution is mu = shape1/(shape1+shape2), which therefore is the mean or the probability of success.

See betabinomial and betabinomialff, the VGAM family functions for estimating the parameters, for the formula of the probability density function and other details.

For the inflated beta-binomial distribution, the probability mass function is $$P(Y = y) = (1 - pstr0 - pstrsize) \times BB(y) + pstr0 \times I[y = 0] + pstrsize \times I[y = size]$$

where $BB(y)$ is the probability mass function of the beta-binomial distribution with the same shape parameters (pbetabinom.ab), pstr0 is the inflated probability at 0 and pstrsize is the inflated probability at 1. The default values of pstr0 and pstrsize mean that these functions behave like the ordinary Betabinom when only the essential arguments are inputted.

Note

pzoibetabinom, pzoibetabinom.ab, pbetabinom and pbetabinom.ab can be particularly slow. The functions here ending in .ab are called from those functions which don't. The simple transformations $\mu=\alpha / (\alpha + \beta)$ and $\rho=1/(1 + \alpha + \beta)$ are used, where $\alpha$ and $\beta$ are the two shape parameters.

Warning

Setting rho = 1 is not recommended, however the code may be modified in the future to handle this special case.

Examples

set.seed(1); rbetabinom(10, 100, prob = 0.5)
#>  [1] 52 46 60 49 52 51 63 52 47 38
set.seed(1);     rbinom(10, 100, prob = 0.5)  # The same as rho = 0
#>  [1] 52 46 60 49 52 51 63 52 47 38

if (FALSE)  N <- 9; xx <- 0:N; s1 <- 2; s2 <- 3
#> Error: object 'N' not found
dy <- dbetabinom.ab(xx, size = N, shape1 = s1, shape2 = s2)
#> Error: object 'xx' not found
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)
#> Error: object 'dy' not found
sum(dy * xx)  # Check expected values are equal
#> Error: object 'dy' not found
sum(dbinom(xx, size = N, prob = s1 / (s1+s2)) * xx)
#> Error: object 'xx' not found
# Should be all 0:
cumsum(dy) - pbetabinom.ab(xx, N, shape1 = s1, shape2 = s2)
#> Error: object 'dy' not found

y <- rbetabinom.ab(n = 1e4, size = N, shape1 = s1, shape2 = s2)
#> Error: object 'N' not found
ty <- table(y)
#> Error: object 'y' not found
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))
#> Error: object 'dy' not found

N <- 1e5; size <- 20; pstr0 <- 0.2; pstrsize <- 0.2
kk <- rzoibetabinom.ab(N, size, s1, s2, pstr0, pstrsize)
#> Error: object 's1' not found
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))
#> Error: object 'kk' not found
sum(kk == 0) / N  # Proportion of 0
#> Error: object 'kk' not found
sum(kk == size) / N  # Proportion of size
#> Error: object 'kk' not found
lines(0 : size,
      dbetabinom.ab(0 : size, size, s1, s2), col = "orange")
#> Error: object 's1' not found
lines(0 : size, col = "green", type = "b",
      dzoibetabinom.ab(0 : size, size, s1, s2, pstr0, pstrsize))
#> Error: object 's1' not found
 # \dontrun{}