# Visualize bivariate binomial distribution

Question: what does a bivariate binomial distribution look like in 3-dimensional space?

Below is the specific function that I would like to visualize for various values of the parameters; namely, $n$, $p_{1}$, and $p_{2}$.

$$f(x_{1},x_{2}) = \frac{n!}{x_{1}!x_{2}!}p_{1}^{x_{1}}p_{2}^{x_{2}}, \qquad x_{1}+x_{2}=n, \quad p_{1}+p_{2}=1.$$

Notice that there are two constraints; $x_{1}+x_{2}=n$ and $p_{1}+p_{2}=1$. In addition, $n$ is a positive integer, say, $5$.

In have made two attempts to plot the function using LaTeX (TikZ/PGFPLOTS). In doing so, I get the graphs below for the following values: $n=5$, $p_{1}=0.1$ and $p_{2}=0.9$, and, $n=5$, $p_{1}=0.4$ and $p_{2}=0.6$, respectively. I haven't been successful at implementing the constraint on the domain values; $x_{1}+x_{2}=n$, so I'm a bit stumped.

A visualization produced in any language would do fine (R, MATLAB, etc.), but I'm working in LaTeX with TikZ/PGFPLOTS.

First Attempt

$n=5$, $p_{1}=0.1$ and $p_{2}=0.9$ Second Attempt

$n=5$, $p_{1}=0.4$ and $p_{2}=0.6$ Edit:

For reference, here is an article containing some graphs. Title of paper is "A new bivariate binomial distribution" by Atanu Biswasa and Jing-Shiang Hwang. Statistics & Probability Letters 60 (2002) 231–240.

Edit 2: For clarity, and in response to @GlenB in the comments, below is a snapshot of how the distribution has been presented to me in my book. The book does not refer to degenerate / non-degenerate cases and so on. It simply presents it like that and I sought to visualize it. Cheers! Also, as pointed out by @JohnK, there is likely to be a typo with regard to x1+x1=1, which he suggests should be x1+x1=n. Image of equation from:

Spanos, A (1986) Statistical foundations of econometric modelling. Cambridge University Press

• But it shouldn't be a continuous, should it? Both random variables are discrete. – JohnK Jan 7 '16 at 19:32
• So x1 & x2 are independent, is that right? You need a pseudo-3D plot? Would a heatmap be acceptable? – gung - Reinstate Monica Jan 7 '16 at 19:32
• something like this? – Antoni Parellada Jan 7 '16 at 19:38
• @JohnK If $x_1+x_2=n$ and $p_1+p_2=1$ you're dealing with $X_1\sim \text{Binomial}(n,p_1)$ (and $X_2$ is simply $n-X_1$). This is univariate binomial (or, considered as bivariate, it's degenerate). – Glen_b -Reinstate Monica Jan 8 '16 at 7:53
• You don't have a specification for a bivariate binomial in your question. (There's more than one way to specify a bivariate distribution that could plausibly be called "binomial". You don't have any of them, though your degenerate one would be a special case of some of them.) ... the drawings in your Biswasa&Hwang reference are not suitable displays of a discrete bivariate pmf. In short, your question lacks anything to draw, and your reference is useful mainly as an example of what to avoid. – Glen_b -Reinstate Monica Jan 8 '16 at 9:48

There are two pieces to this: first you need to figure out what the individual probabilities are, then you need to plot them somehow.

A binomial PMF is just a set of probabilities over a number of 'successes'. A bivariate binomial PMF will be a set of probabilities over a grid of possible combinations of 'successes'. In your case, you have $n_i = n_j = 5$, so (bearing in mind that $0$ successes is a possibility) there are $6\times 6 = 36$ possible outcomes in the grid / bivariate binomial distribution.

We can first calculate the marginal binomial PMFs, because that is so straightforward. Since the variables are independent, each joint probability will just be the product of the marginal probabilities; this is matrix algebra. Here I demonstrate this process using R code:

b1 = dbinom(0:5, size=5, prob=0.1);  sum(b1)  #  1
b9 = dbinom(0:5, size=5, prob=0.9);  sum(b9)  #  1
b4 = dbinom(0:5, size=5, prob=0.4);  sum(b4)  #  1
b6 = dbinom(0:5, size=5, prob=0.6);  sum(b6)  #  1

b19 = b1%o%b9;  sum(b19)  #  1
rownames(b19) <- colnames(b19) <- as.character(0:5)
round(b19, 6)
#       0        1        2        3        4        5
# 0 6e-06 0.000266 0.004783 0.043047 0.193710 0.348678
# 1 3e-06 0.000148 0.002657 0.023915 0.107617 0.193710
# 2 1e-06 0.000033 0.000590 0.005314 0.023915 0.043047
# 3 0e+00 0.000004 0.000066 0.000590 0.002657 0.004783
# 4 0e+00 0.000000 0.000004 0.000033 0.000148 0.000266
# 5 0e+00 0.000000 0.000000 0.000001 0.000003 0.000006
b46 = b4%o%b6;  sum(b46)  #  1
rownames(b46) <- colnames(b46) <- as.character(0:5)
round(b46, 3)
#       0     1     2     3     4     5
# 0 0.001 0.006 0.018 0.027 0.020 0.006
# 1 0.003 0.020 0.060 0.090 0.067 0.020
# 2 0.004 0.027 0.080 0.119 0.090 0.027
# 3 0.002 0.018 0.053 0.080 0.060 0.018
# 4 0.001 0.006 0.018 0.027 0.020 0.006
# 5 0.000 0.001 0.002 0.004 0.003 0.001


At this point, we have the two requisite matrices of probabilities. We just need to decide how we want to plot them. To be honest, I am not a big fan of 3D bars charts. Because R seems to agree with me, I made these plots in Excel:

b19: b46: • Thank you for presentation plus R code. This leads me to ask about x1+x2=n. If this condition holds, should there only be a single line of pillars as presented here: reference.wolfram.com/language/ref/MultinomialDistribution.html The wolfram graph I assume is what @Glen_b has referred to as degenerate case? Does this mean that you've presented the non-degenerate case? – Graeme Walsh Jan 8 '16 at 9:48
• GraemeWalsh, my presentation does not show a bivariate binomial where x1+x2=n. As @Glen_b discussed extensively in the comments & his answer, I wouldn't really call that a "bivariate binomial distribution" w/o qualifying it. Moreover, it would mean that x1 & x2 are not independent, as you said in your response comment, but perfectly dependent. In truth, I didn't notice that this was such a bizarre variant (you can blame me for not reading closely enough). As Glen_b showed, that version would be a single line of pillars. What I presented was the non-degenerate case. – gung - Reinstate Monica Jan 8 '16 at 16:24
• @gung I like your new plots. I think your discussion covers the degenerate case just fine ("you need to figure out what the individual probabilities are" really says everything; the actual calculations for the degenerate case are trivial); I just carried out those trivial calculations. – Glen_b -Reinstate Monica Jan 9 '16 at 8:17

gung's answer is a good answer for an actual bivariate binomial, explaining the issues well (I'd recommend accepting it as a good answer to the title question, most likely to be useful to others).

The mathematical object you actually present in your edit is really a univariate scaled binomial. Here $x_1$ is not the value taken by the binomial count but by the proportion (the binomial divided by $n$).

So let's define things properly. Note that no definition of the random variable is actually offered, so we're left with some guesswork.

Let $Y_1\sim \text{binomial}(n,p_1),\:$ Note that when we give a mathematical formula for $P(Y_1=y_1)$ it's necessary what values $y_1$ can take, so $y_1=0,1,...,n$. Let $X_1=Y_1/n$, and note that $x_1=0,\frac16,\frac26,...,1$.

Then the equation you give is the pmf for $P(X_1=x_1)$ (noting that $x_2=n-x_1$ and $p_2=1-p_1$).

For $n=6,p_1=0.3$, it looks like this: We can put $x_2$ values on the above plot quite readily, simply by putting a second set of labels under the $x_1$ values equal to $1-x_1$ (perhaps in a different colour) to indicate the value taken by $x_2$.

We could regard it as a (scaled) degenerate bivariate binomial: but it's a bit of a stretch to really call what's defined in the book a bivariate binomial, (since it's effectively a univariate binomial).

On the assumption that someone will want to generate a similar plot to the 3D one, this little bit of (R) code gets quite close to the second plot above:

y = 0:6
x1 = y/6
x2 = 1-x1
p = dbinom(y,6,.3)
scatterplot3d(x1,x2,p,grid=TRUE, box=FALSE, cex.lab=1.2,
color=3, cex.main=1.4,pch=21,bg=1,, type="h",angle=120,
main="degenerate scaled binomial", ylab="x2", xlab="x1",
zlab="prob")


(You need the scatterplot3d package which contains the function of the same name.)

A "true" (non-degenerate) bivariate binomial has variation in both variables at once. Here's an example of one particular kind of bivariate binomial (not independent in this case). I resorted to using different colours in the plot because it's too easy to get lost in the forest of "sticks" otherwise. There are many ways to get an object that you might call a bivariate binomial; this particular kind is one where you have $X\sim\text{bin}(n_0,p)$,$Y\sim\text{bin}(n_y,p)$,$Z\sim\text{bin}(n_z,p)$ (all independent), then let $X_1=X+Y$ and $X_2=X+Z$.

This yields binomial $X_1$ and $X_2$ which are correlated (but has the disadvantage that it doesn't produce negative correlations).

An expression for the pmf of this particular kind of bivariate binomial distribution is given in Hamdan, 1972  but I didn't use that calculation; one can easily do direct computation (numeric convolution). In this particular case $n_0$ was 4 and $n_y$ and $n_z$ were only 2 each so direct numeric computation across the whole grid (49 values in the final result) is not difficult or onerous. You start with a degenerate bivariate (both dimensions $=X$) similar to the degenerate one pictured above (but smaller and on the "main diagonal" - $x_1=x_2$ rather than the antidiagonal ($x_1+x_2=n$) and then add the independent components, spreading the probability along and out from the diagonal.

: Hamdan, M.A. (1972),
"Canonical Expansion of the Bivariate Binomial Distribution with Unequal Marginal Indices"
International Statistical Review, 40:3 (Dec.), pp. 277-280

• Nice. It is also worth noting that in this case $corr(X_1, X_2) = -1$ – JohnK Jan 8 '16 at 11:46
• Glen_b. Thank you very much. Pointing out that the math object I presented (that was presented to me!) is a (scaled) degenerate bivariate binomial has been very helpful! I did not know this from the beginning. Lastly, an elementary request! Would it be possible for you to be explicit (by way of math notation) about how you define a true or actual bivariate binomial? That would be useful, I think. – Graeme Walsh Jan 8 '16 at 15:07
• @Graeme As I mention already in comments (/answer), there are many ways to get an object that you might call a bivariate binomial (indeed the title of the Biswasa and Hwang reference in your question tells you as much). This is not unique to the binomial of course, there are many available bivariate generalizations of a lot of the more commonly used univariate distributions. The "particular kind of bivariate binomial" I gave in my answer is one where you have $X\sim\text{bin}(n_0,p)$,$Y\sim\text{bin}(n_y,p)$,$Z\sim\text{bin}(n_z,p)$ (all independent), then let $X_1=X+Y$ and $X_2=X+Z$. ... ctd – Glen_b -Reinstate Monica Jan 8 '16 at 20:25
• ctd ...This yields binomial $X_1$ and $X_2$ which are correlated but has the disadvantage that it doesn't produce negative correlations, so it's not as useful as some other formulations of bivariate binomial for general bivariate modelling. Typically when you generalize a univariate distribution family to a bivariate one, you have to choose which properties you want most and which you can afford to give up, and those choices will lead to different choices of bivariate families. [The normal distribution is unusual - there's one "obvious" generalization with pretty much everything we'd want.] – Glen_b -Reinstate Monica Jan 8 '16 at 20:33
• @Graeme ... I plan to add some more details. – Glen_b -Reinstate Monica Jan 9 '16 at 2:10

Mathematica is now quite strong in such things - it has the solution of your problem right in documentation. With little additions I've made a model to play around (with p = p1 = 0.4 for better visual presentation). That is how interface looks and how it can be controlled. Snippet

Manipulate[
Grid[{
{DiscretePlot3D[
PDF[MultinomialDistribution[n, {p, 1 - p}], {x, y}], {x, 0,
n}, {y, 0, n}, PlotLabel -> Row[{"n = ", n}],
ExtentSize -> Right],

DiscretePlot3D[
CDF[MultinomialDistribution[n, {p, 1 - p}], {x, y}], {x, 0,
n}, {y, 0, n}, PlotLabel -> Row[{"n = ", n}],
ExtentSize -> Right]}
}]
,
{{n, 5}, 1, 20, 1, Appearance -> "Labeled"},
{{p, 0.4}, 0.1, 0.9},
TrackedSymbols -> True
]


The main thing here is PDF[MultinomialDistribution[n, {p, 1 - p}], {x, y}], which is selfexplanatory, I think. Multinomial just mean that you may take a lot of distributions with each pi for respective variable. The simple form is BinomialDistribution. Of course, I could make it manually, but the rule is if you have a build-in function - you should use it.