BetaBinom function

Beta-binomial distribution

Beta-binomial distribution

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

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[Xx]P[X \le x]

    otherwise, P[X>x]P[X > x].

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

Details

If p Beta(α,β)p ~ Beta(\alpha, \beta) and X Binomial(n,p)X ~ Binomial(n, p), then X BetaBinomial(n,α,β)X ~ BetaBinomial(n, \alpha, \beta).

Probability mass function

f(x)=(nx)B(x+α,nx+β)B(α,β)f(x)=choose(n,x)B(x+α,nx+β)/B(α,β) f(x) = {n \choose x} \frac{\mathrm{B}(x+\alpha, n-x+\beta)}{\mathrm{B}(\alpha, \beta)}f(x) = choose(n, x) * B(x+\alpha, n-x+\beta) / B(\alpha, \beta)

Cumulative distribution function is calculated using recursive algorithm that employs the fact that Γ(x)=(x1)!\Gamma(x) = (x - 1)!, and c("\n\n", "B(x,y)=(Gamma(x)Gamma(y))/Gamma(x+y)\nB(x, y) = (\\Gamma(x)\\Gamma(y))/\\Gamma(x+y)\n"), and that c("\n\n", "choose(n,k)=prod((n+1(1:k))/(1:k))\nchoose(n,k) = prod((n+1-(1:k))/(1:k))\n"). This enables re-writing probability mass function as

f(x)=(i=1xn+1ii)(α+x1)!(β+nx1)!(α+β+n1)!B(α,β)f(x)=prod((n+1(1:x))/(1:x))(((α+x1)!(β+nx1)!)/((α+β+n+1)!))/B(α,β) 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)}f(x) = prod((n+1-(1:x))/(1:x)) * (((\alpha+x-1)!*(\beta+n-x-1)!)/((\alpha+\beta+n+1)!)) / B(\alpha, \beta)

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

f(x+1)=(i=1xn+1ii)n+1x+1x+1(α+x1)!(α+x)(β+nx1)!(β+nx)1(α+β+n1)!(α+β+n)B(α,β)f(x+1)=prod((n+1(1:x))/(1:x))((n+1x+1)/(x+1))(((α+x1)!(α+x)(β+nx1)!/(β+nx))/((α+β+n+1)!(α+β+n)))/B(α,β) 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)}f(x+1) = prod((n+1-(1:x))/(1:x)) * ((n+1-x+1)/(x+1)) * (((\alpha+x-1)!*(\alpha+x)*(\beta+n-x-1)!/(\beta+n-x))/((\alpha+\beta+n+1)!*(\alpha+\beta+n))) / B(\alpha, \beta)

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

F(x)=k=0xf(k)F(x)=f(0)+...+f(x) F(x) = \sum_{k=0}^x f(k)F(x) = f(0)+...+f(x)

Examples

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)

See Also

Beta, Binomial

  • Maintainer: Tymoteusz Wolodzko
  • License: GPL-2
  • Last published: 2023-11-30