# BetaBinom: Beta-binomial distribution

## Description

Probability mass function and random generation for the beta-binomial distribution.

## Usage

dbbinom(x, size, alpha = 1, beta = 1, log = FALSE)pbbinom(q, size, alpha = 1, beta = 1, lower.tail = TRUE, log.p = FALSE)rbbinom(n, size, alpha = 1, beta = 1)

## Arguments

x, q

vector of quantiles.

size

number of trials (zero or more).

alpha, beta

non-negative parameters of the beta distribution.

log, log.p

logical; if TRUE, probabilities p are given as log(p).

lower.tail

logical; if TRUE (default), probabilities are $$P[X \le x]$$ otherwise, $$P[X > x]$$.

n

number of observations. If length(n) > 1, the length is taken to be the number required.

## Details

If $$p \sim \mathrm{Beta}(\alpha, \beta)$$ and $$X \sim \mathrm{Binomial}(n, p)$$, then $$X \sim \mathrm{BetaBinomial}(n, \alpha, \beta)$$.

Probability mass function $$f(x) = {n \choose x} \frac{\mathrm{B}(x+\alpha, n-x+\beta)}{\mathrm{B}(\alpha, \beta)}$$

Cumulative distribution function is calculated using recursive algorithm that employs the fact that $$\Gamma(x) = (x - 1)!$$, and $$\mathrm{B}(x, y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}$$, and that $${n \choose k} = \prod_{i=1}^k \frac{n+1-i}{i}$$. This enables re-writing probability mass function as

$$f(x) = \left( \prod_{i=1}^x \frac{n+1-i}{i} \right) \frac{\frac{(\alpha+x-1)!\,(\beta+n-x-1)!}{(\alpha+\beta+n-1)!}}{\mathrm{B}(\alpha,\beta)}$$

what makes recursive updating from $$x$$ to $$x+1$$ easy using the properties of factorials

$$f(x+1) = \left( \prod_{i=1}^x \frac{n+1-i}{i} \right) \frac{n+1-x+1}{x+1} \frac{\frac{(\alpha+x-1)! \,(\alpha+x)\,(\beta+n-x-1)! \, (\beta+n-x)^{-1}}{(\alpha+\beta+n-1)!\,(\alpha+\beta+n)}}{\mathrm{B}(\alpha,\beta)}$$

and let's us efficiently calculate cumulative distribution function as a sum of probability mass functions

$$F(x) = \sum_{k=0}^x f(k)$$

Beta, Binomial

## Examples

# NOT RUN {
x <- rbbinom(1e5, 1000, 5, 13)
xx <- 0:1000
hist(x, 100, freq = FALSE)
lines(xx-0.5, dbbinom(xx, 1000, 5, 13), col = "red")
hist(pbbinom(x, 1000, 5, 13))
xx <- seq(0, 1000, by = 0.1)
plot(ecdf(x))
lines(xx, pbbinom(xx, 1000, 5, 13), col = "red", lwd = 2)

# }