Betabinom: The Beta-Binomial Distribution

Description

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

Usage

dbetabinom(x, size, prob, rho = 0, log = FALSE)
pbetabinom(q, size, prob, rho, log.p = FALSE)
rbetabinom(n, size, prob, rho = 0)
dbetabinom.ab(x, size, shape1, shape2, log = FALSE, .dontuse.prob = NULL)
pbetabinom.ab(q, size, shape1, shape2, log.p = FALSE)
rbetabinom.ab(n, size, shape1, shape2, .dontuse.prob = NULL)

Arguments

x, q

vector of quantiles.

size

number of trials.

number of observations. Must be a positive integer of length 1.

prob

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

rho

the correlation parameter $\rho$. Usually must be in the unit open interval $(0,1)$, however, the value 0 is sometimes supported (if so then it corresponds to the usual binomial distribution).

shape1, shape2

the two (positive) shape parameters of the standard beta distribution. They are called a and b in beta respectively.

log, log.p

Logical. If TRUE then all probabilities p are given as log(p).

.dontuse.prob

An argument that should be ignored and unused.

Value

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

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 betabinomial.ab, the VGAM family functions for estimating the parameters, for the formula of the probability density function and other details.

Examples

Run this code

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)

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","red"), 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),
        ") (red)", sep = ""), cex.main = 0.8,
        names.arg = as.character(xx))