# The Beta-Binomial Distribution Density, distribution function, and random generation for the beta-binomial distribution and the inflated beta-binomial distribution. ```r 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. ## Returns `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(s) 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]%F(Y = y) =(1 - pstr0 - pstrsize) * BB(y) + pstr0 * I[y = 0] +pstrsize * 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. ## See Also `Extbetabinom`, `betabinomial`, `betabinomialff`, `Zoabeta`, `Beta`. ## Examples ```r set.seed(1); rbetabinom(10, 100, prob = 0.5) set.seed(1); rbinom(10, 100, prob = 0.5) # The same as 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) # Should be all 0: cumsum(dy) - pbetabinom.ab(xx, N, shape1 = s1, shape2 = s2) 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)) ## End(Not run) ```

betabinomUC function