Poisson binomial distribution

The probability of having k successful trials out of a total of n can be written as the sum

{{Cite journal

| volume = 3

| issue = 2

| pages = 295–312

| last = Wang

| first = Y. H.

| title = On the number of successes in independent trials

| journal = Statistica Sinica

| year = 1993

| url = http://www3.stat.sinica.edu.tw/statistica/oldpdf/A3n23.pdf

}}

:$\backslash Pr(K=k)\; =\; \backslash sum\backslash limits\_\{A\backslash in\; F\_k\}\; \backslash prod\backslash limits\_\{i\backslash in\; A\}\; p\_i\; \backslash prod\backslash limits\_\{j\backslash in\; A^c\}\; (1-p\_j)$

where $F\_k$ is the set of all subsets of k integers that can be selected from $\backslash \{1,2,3,...,n\backslash \}$. For example, if n = 3, then $F\_2=\backslash left\backslash \{\; \backslash \{1,2\backslash \},\backslash \{1,3\backslash \},\backslash \{2,3\backslash \}\; \backslash right\backslash \}$. $A^c$ is the complement of $A$, i.e. $A^c\; =\backslash \{1,2,3,\backslash dots,n\backslash \}\backslash smallsetminus\; A$.

$F\_k$ will contain $n!/((n-k)!k!)$ elements, the sum over which is infeasible to compute in practice unless the number of trials n is small (e.g. if n = 30, $F\_\{15\}$ contains over 10²⁰ elements). However, there are other, more efficient ways to calculate $\backslash Pr(K=k)$.

As long as none of the success probabilities are equal to one, one can calculate the probability of k successes using the recursive formula

{{Cite journal

| volume = 27

| issue = 3

| pages = 123–124

| last = Shah

| first = B. K.

| title = On the distribution of the sum of independent integer valued random variables

| journal = American Statistician

| year = 1994

| jstor=2683639

}}

{{Cite journal

| volume = 81

| issue = 3

| pages = 457

| last = Chen

| first = X. H.

|author2=A. P. Dempster |author3=J. S. Liu

| title = Weighted finite population sampling to maximize entropy

| journal = Biometrika

| year = 1994

| url = http://www.people.fas.harvard.edu/~junliu/TechRept/94folder/cdl94.pdf

| doi=10.1093/biomet/81.3.457

}}

:$\backslash Pr\; (K=k)=\; \backslash begin\{cases\}$

\prod\limits_{i=1}^n (1-p_i) & k=0 \\

\frac{1}{k} \sum\limits_{i=1}^k (-1)^{i-1}\Pr (K=k-i)T(i) & k>0 \\

\end{cases}

where

: $T(i)=\backslash sum\backslash limits\_\{j=1\}^n\; \backslash left(\; \backslash frac\{p\_j\}\{1-p\_j\}\; \backslash right)^i.$

The recursive formula is not numerically stable, and should be avoided if $n$ is greater than approximately 20.

An alternative is to use a divide-and-conquer algorithm: if we assume $n\; =\; 2^b$ is a power of two, denoting by $f(p\_\{i:j\})$ the Poisson binomial of $p\_i,\; \backslash dots,\; p\_j$ and $*$ the convolution operator, we have $f(p\_\{1:2^b\})\; =\; f(p\_\{1:2^\{b-1\}\})*f(p\_\{2^\{b-1\}+1:2^b\})$.

More generally, the probability mass function of a Poisson binomial can be expressed as the convolution of the vectors $P\_1,\; \backslash dots,\; P\_n$ where $P\_i=[1-p\_i,p\_i]$. This observation leads to the direct convolution (DC) algorithm for computing $\backslash Pr\; (K=0)$ through $\backslash Pr\; (K=n)$:

// PMF and nextPMF begin at index 0

function DC($p\_1,\; \backslash dots,\; p\_n$) is

declare new PMF array of size 1

PMF[0] = [1]

for i = 1 to $n$ do

declare new nextPMF array of size i + 1

nextPMF[0] = (1 - $p\_i$) * PMF[0]

nextPMF[i] = $p\_i$ * PMF[i - 1]

for k = 1 to i - 1 do

nextPMF[k] = $p\_i$ * PMF[k - 1] + (1 - $p\_i$) * PMF[k]

repeat

PMF = nextPMF

repeat

return PMF

end function

$\backslash Pr\; (K=k)$will be found in PMF[k]. DC is numerically stable, exact, and, when implemented as a software routine, exceptionally fast for $n\; \backslash leq\; 2000$. It can also be quite fast for larger $n$, depending on the distribution of the $p\_i$.

Another possibility is using the discrete Fourier transform.

{{Cite journal

| volume = 46

| issue = 2

| pages = 803–817

| last = Fernandez

| first = M.

|author2=S. Williams

| title = Closed-Form Expression for the Poisson-Binomial Probability Density Function

| journal = IEEE Transactions on Aerospace and Electronic Systems

| year = 2010

| doi = 10.1109/TAES.2010.5461658

| bibcode = 2010ITAES..46..803F

| s2cid = 1456258

}}

:$\backslash Pr\; (K=k)=\backslash frac\{1\}\{n+1\}\; \backslash sum\_\{\backslash ell=0\}^n\; C^\{-lk\}\; \backslash prod\_\{m=1\}^n\; \backslash left(\; 1+(C^\backslash ell-1)\; p\_m\; \backslash right)$

where $C=\backslash exp\; \backslash left(\; \backslash frac\{2i\backslash pi\; \}\{n+1\}\; \backslash right)$ and $i=\backslash sqrt\{-1\}$.

Still other methods are described in "Statistical Applications of the Poisson-Binomial and conditional Bernoulli distributions" by Chen and Liu{{Cite journal

| volume = 7

| pages = 875–892

| last = Chen

| first = S. X.

|author2=J. S. Liu

| title = Statistical Applications of the Poisson–Binomial and conditional Bernoulli distributions

| journal = Statistica Sinica

| year = 1997

| url = http://www3.stat.sinica.edu.tw/statistica/password.asp?vol=7&num=4&art=4

}} and in "A simple and fast method for computing the Poisson binomial distribution function" by Biscarri et al.{{Cite journal |last1=Biscarri |first1=William |last2=Zhao |first2=Sihai Dave |last3=Brunner |first3=Robert J. |date=2018-06-01 |title=A simple and fast method for computing the Poisson binomial distribution function |url=https://www.sciencedirect.com/science/article/pii/S0167947318300082 |journal=Computational Statistics & Data Analysis |volume=122 |pages=92–100 |doi=10.1016/j.csda.2018.01.007 |issn=0167-9473}}

The cumulative distribution function (CDF) can be expressed as:

: $\backslash Pr(K\; \backslash leq\; k)\; =\; \backslash sum^\{k\}\_\{\backslash ell=0\}\; \backslash sum\_\{A\backslash in\; F\_\backslash ell\}\; \backslash prod\_\{i\backslash in\; A\}\; p\_i\; \backslash prod\_\{j\backslash in\; A^c\}\; (1-p\_j),$

where $F\_\backslash ell$ is the set of all subsets of size $\backslash ell$ that can be selected from $\backslash \{1,2,3,\backslash ldots,n\backslash \}$.

It can be computed by invoking the DC function above, and then adding elements $0$ through $k$ of the returned PMF array.

Since a Poisson binomial distributed variable is a sum of n independent Bernoulli distributed variables, its mean and variance will simply be sums of the mean and variance of the n Bernoulli distributions:

: $\backslash mu\; =\; \backslash sum\backslash limits\_\{i=1\}^n\; p\_i$

: $\backslash sigma^2\; =\backslash sum\backslash limits\_\{i=1\}^n\; (1-p\_i)\; p\_i$

There is no simple formula for the entropy of a Poisson binomial distribution, but the entropy is bounded above by the entropy of a binomial distribution with the same number parameter and the same mean. Therefore, the entropy is also bounded above by the entropy of a Poisson distribution with the same mean.{{Cite journal

| volume = 47 | issue = 5

| pages = 2039–2041

| last = Harremoës

| first = P.

| title = Binomial and Poisson distributions as maximum entropy distributions

| journal = IEEE Transactions on Information Theory

| year = 2001

| url = https://web.archive.org/web/20160116154755id_/http://yaroslavvb.com/papers/harremoes-binomial.pdf

| doi=10.1109/18.930936

}}

The Shepp–Olkin concavity conjecture, due to Lawrence Shepp and Ingram Olkin in 1981, states that the entropy of a Poisson binomial distribution is a concave function of the success probabilities $p\_1,\; p\_2,\; \backslash dots\; ,\; p\_n$.{{Cite encyclopedia

|editor1-first = J.

|editor1-last = Gani

|editor2-first = V.K.

|editor2-last = Rohatgi

|author1-last = Shepp

|author1-first = Lawrence

|author2-last = Olkin

|author2-first = Ingram

|title = Entropy of the sum of independent Bernoulli random variables and of the multinomial distribution

|encyclopedia = Contributions to probability: A collection of papers dedicated to Eugene Lukacs

|pages = 201–206

|publisher = Academic Press

|location = New York

|year = 1981

|mr = 0618689

|url = https://books.google.com/books?id=9qziBQAAQBAJ

|isbn = 0-12-274460-8

}} This conjecture was proved by Erwan Hillion and Oliver Johnson in 2015.{{Cite journal |last1=Hillion|first1=Erwan|last2=Johnson|first2=Oliver|date=2015-03-05|title=A proof of the Shepp–Olkin entropy concavity conjecture |journal=Bernoulli|volume = 23|issue=4B|pages = 3638–3649 | arxiv=1503.01570 |doi=10.3150/16-BEJ860|s2cid=8358662}} The Shepp–Olkin monotonicity conjecture, also from the same 1981 paper, is that the entropy is monotone increasing in $p\_i$, if all $p\_i\; \backslash leq\; 1/2$. This conjecture was also proved by Hillion and Johnson, in 2019.{{Cite journal |last1=Hillion|first1=Erwan|last2=Johnson|first2=Oliver|date=2019-11-09|title=A proof of the Shepp–Olkin entropy monotonicity conjecture |journal=Electronic Journal of Probability| volume = 24 |number=126 | pages = 1–14 |doi=10.1214/19-EJP380|doi-access=free|arxiv=1810.09791}}

The probability that a Poisson binomial distribution gets large, can be bounded using its moment generating function as follows (valid when $s\; \backslash geq\; \backslash mu$ and for any $t>0$):

: $$

\begin{align}

\Pr[S\ge s]

&\le \exp(-st)\operatorname E\left[\exp\left[t \sum_i X_i\right]\right]

\\&= \exp(-st)\prod_i (1-p_i+e^t p_i)

\\&= \exp\left(-st + \sum_i \log\left( p_i (e^t - 1) + 1\right)\right)

\\&\le \exp\left(-st + \sum_i \log\left( \exp(p_i(e^t-1))\right)\right)

\\&= \exp\left(-st + \sum_i p_i(e^t-1)\right)

\\&= \exp\left(s-\mu-s\log\frac{s}{\mu}\right),

\end{align}

where we took $t=\backslash log\backslash left(s/\backslash mu\backslash right)$. This is similar to the tail bounds of a binomial distribution.

A Poisson binomial distribution $PB$ can be approximated by a binomial distribution $B$ where $\backslash mu$, the mean of the $p\_i$, is the success probability of $B$. The variances of $PB$ and $B$ are related by the formula

: $\backslash operatorname\{Var\}(PB)=\backslash operatorname\{Var\}(B)-\; \backslash sum\_\{i=1\}^n\; (p\_i-\backslash mu)^2$

As can be seen, the closer the $p\_i$ are to $\backslash mu$, that is, the more the $p\_i$ tend to homogeneity, the larger $PB$'s variance. When all the $p\_i$are equal to $\backslash mu$, $PB$ becomes $B$, $\backslash operatorname\{Var\}(PB)=\backslash operatorname\{Var\}(B)$, and the variance is at its maximum.

Ehm has determined bounds for the total variation distance of $PB$ and $B$, in effect providing bounds on the error introduced when approximating $PB$ with $B$. Let $\backslash nu\; =\; 1\; -\; \backslash mu$ and $d(PB,B)$ be the total variation distance of $PB$ and $B$. Then

: $d(PB,B)\backslash le(1-\backslash mu^\{n+1\}-\backslash nu^\{n+1\})\; \backslash frac\{\backslash sum\_\{i=1\}^n\; (p\_i-\backslash mu)^2\}\{((n+1)\backslash mu\backslash nu)\}$

: $d(PB,B)\backslash ge\; C\backslash min\backslash left\backslash \{\backslash ,1,\backslash frac\{1\}\{n\backslash mu\backslash nu\}\backslash ,\backslash right\backslash \}\; \backslash sum\_\{i=1\}^n\; (p\_i-\backslash mu)^2$

where $C\backslash ge\backslash frac\{1\}\{124\}$.

$d(PB,B)$ tends to 0 if and only if $\backslash operatorname\{Var\}(PB)/\backslash operatorname\{Var\}(B)$ tends to 1.{{Cite journal |last=Ehm |first=Werner |date=1991-01-01 |title=Binomial approximation to the Poisson binomial distribution |url=https://dx.doi.org/10.1016/0167-7152%2891%2990170-V |journal=Statistics & Probability Letters |volume=11 |issue=1 |pages=7–16 |doi=10.1016/0167-7152(91)90170-V |issn=0167-7152|url-access=subscription }}

A Poisson binomial distribution $PB$ can also be approximated by a Poisson distribution $Po$ with mean $\backslash lambda=\backslash sum\_\{i=1\}^n\; p\_i$

. Barbour and Hall have shown that

: $\backslash frac\{1\}\{32\}\backslash min\backslash left\backslash \{\backslash ,\backslash frac\{1\}\{\backslash lambda\},1\backslash ,\backslash right\backslash \}\; \backslash sum\_\{i=1\}^n\; p\_i^2\backslash le\; d(PB,\; Po)\backslash le\backslash frac\{1-\backslash varepsilon^\{-\backslash lambda\}\}\; \backslash lambda\; \backslash sum\_\{i=1\}^n\; p\_i^2$

where $d(PB,B)$ is the total variation distance of $PB$ and $Po$.{{Cite journal

| last1=Barbour | first1=A. D.

| last2=Hall | first2=Peter

| date=1984

| title=On the Rate of Poisson Convergence | url=https://www.zora.uzh.ch/id/eprint/23054/11/Barbour_1984V.pdf

| journal=Mathematical Proceedings of the Cambridge Philosophical Society

| volume=95

| issue=3

| pages=473–480

| doi=10.1017/S0305004100061806| doi-broken-date=24 December 2024

}} It can be seen that the smaller the $p\_i$, the better $Po$ approximates $PB$.

As $\backslash operatorname\{Var\}(Po)=\backslash lambda=\backslash sum\_\{i=1\}^n\; p\_i$

and $\backslash operatorname\{Var\}(PB)\; =\backslash sum\backslash limits\_\{i=1\}^n\; p\_i-\backslash sum\backslash limits\_\{i=1\}^n\; p\_i^2$, $\backslash operatorname\{Var\}(\backslash mathrm\{Po\})>\backslash operatorname\{Var\}(PB)$

; so a Poisson binomial distribution's variance is bounded above by a Poisson distribution with $\backslash lambda=\backslash sum\_\{i=1\}^n\; p\_i$, and the smaller the $p\_i$, the closer $\backslash operatorname\{Var\}(\backslash mathrm\{Po\})$

will be to $\backslash operatorname\{Var\}(PB)$

.

The reference

{{cite journal |last1=Hong |first1=Yili |title=On computing the distribution function for the Poisson binomial distribution |journal=Computational Statistics & Data Analysis |date=March 2013 |volume=59 |pages=41–51 | doi = 10.1016/j.csda.2012.10.006}}

discusses techniques of evaluating the probability mass function of the Poisson binomial distribution. The following software implementations are based on it:

• An R package [https://CRAN.R-project.org/package=poibin poibin] was provided along with the paper, which is available for the computing of the cdf, pmf, quantile function, and random number generation of the Poisson binomial distribution. For computing the PMF, a DFT algorithm or a recursive algorithm can be specified to compute the exact PMF, and approximation methods using the normal and Poisson distribution can also be specified.
• [https://github.com/tsakim/poibin poibin – Python implementation] – can compute the PMF and CDF, uses the DFT method described in the paper for doing so.

{{Portal|Mathematics}}

• Le Cam's theorem

{{ProbDistributions|discrete-finite}}

Category:Discrete distributions

Category:Factorial and binomial topics