# Examples of variances of random variables¶

Let’s revisit some of the distributions we encountered in the earlier and calculate their variances. We will do it both analytically, and using scipy.stats.

## Example: Variance of a Bernoulli random variable¶

Take a Bernoulli random variable:

$X \sim \text{Bernoulli}(\theta).$

We have already found that this expectation is:

$\mu = \mathbb{E}[X] = \theta.$

To find the variance we are going to use Variance Property 3. For this we need to find the expectation of the square:

$\begin{split} \begin{split} \mathbb{E}[X^2] &= \sum_x x^2 p(x)\\ &= 0^2 p(X=0) + 1^2 p(X=1)\\ &= \theta. \end{split} \end{split}$

So, we have:

$\mathbb{V}[X] = \mathbb{E}[X^2] - \mu^2 = \theta - \theta^2 = \theta (1 - \theta).$

And here is how we can do it using scipy.stats:

import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set(rc={"figure.dpi":100, 'savefig.dpi':300})
sns.set_context('notebook')
sns.set_style("ticks")
from IPython.display import set_matplotlib_formats
set_matplotlib_formats('retina', 'svg')
import numpy as np
import scipy.stats as st

theta = 0.7
X = st.bernoulli(theta)


Now that we have made the random variable we can get its expectation by X.var():

print('V[X] = {0:1.2f}'.format(X.var()))
print('Compare to theta * (1 - theta) = {0:1.2f}'.format(theta * (1 - theta)))

V[X] = 0.21
Compare to theta * (1 - theta) = 0.21


The standard deviation is just the square root of the variance:

$\sigma = \sqrt{\theta (1-\theta)}.$

In scipy.stats you can get it by X.std():

print('std of X = {0:1.2f}'.format(X.std()))

std of X = 0.46


## Example: Variance of a uniform random variable¶

Take

$X \sim U([a,b]).$

Remember that the PDF is:

$p(x) = \frac{1}{b-a},$

when $$x$$ is in $$[a,b]$$ and zero otherwise.

We have already found the expectation and it was given by the mid-point between $$a$$ and $$b$$:

$\mu = \mathbf{E}[X] = \frac{a+b}{2}.$

To find the variance, we first need to find the expectation of the square:

$\begin{split} \begin{split} \mathbb{E}[X^2] &= \int x^2p(x)dx\\ &= \int_a^b x^2 \frac{1}{b-a}dx\\ &= \frac{1}{b-a}\frac{x^3}{3}|_a^b\\ &= \frac{1}{b-a}\frac{b^3 - a^3}{3}\\ &=\frac{b^3 - a^3}{3(b-a)} \end{split} \end{split}$

You can simplify this even more, but we won’t bother. Now, put everything together:

$\begin{split} \begin{split} \mathbf{V}[X] &= \mathbb{E}[X^2] - \mu^2\\ &= \frac{b^3 - a^3}{3(b-a)} - \frac{(a+b)^2}{4}\\ &= \frac{4\cdot (b^3 - a^3)}{4\cdot 3(b-a)} - \frac{3(b-a)\cdot (a+b)^2}{3(b-a)\cdot 4}\\ &= \frac{4\cdot (b^3 - a^3) - 3(b-a)\cdot (a+b)^2}{12(b-a)}\\ &= \frac{4b^3 - 4a^3 - 3(b-a)(a^2 + 2ab+b^2)}{12(b-a)}\\ &= \frac{4b^3 - 4a^3 -3a^2b -6ab^2 - 3b^3 + 3a^3 + 6a^2b + 3ab^2}{12(b-a)}\\ &= \frac{b^3 - a^3 -3a^2b -6ab^2 + 6a^2b + 3ab^2}{12(b-a)}\\ &= \frac{(b-a)^3}{12(b-a)}\\ &= \frac{(b-a)^2}{12}. \end{split} \end{split}$

If you remember your basics mechanics course, this is the second area moment of inertia of a beam about its center of mass. This is not an accident. Mathematically, the variance and the second area moments are exactly the same integrals.

Let’s do it also on scipy.stats:

a = 0
b = 5
X = st.uniform(a, b)
print('V[X] = {0:1.2f}'.format(X.var()))
print('Compare to theoretical answer = {0:1.2f}'.format((b - a) ** 2 / 12))

V[X] = 2.08
Compare to theoretical answer = 2.08


## Example: Variance of a Categorical random variable¶

Take a Categorical random variable:

$X \sim \text{Categorical}(0.1, 0.3, 0.4, 0.2).$

The expectation is:

$\mu = \mathbf{E}[X] = 1.7.$

Again, we are going to invoke Variance Property 3. We need the expectation of the square:

$\begin{split} \begin{split} \mathbf{E}[X^2] &= \sum_x x^2 p(x)\\ &= 0^2\cdot p(X=0) + 1^2 \cdot p(X=1) + 2^2\cdot p(X=2) + 3^2 \cdot p(X=3)\\ &= 0 \cdot 0.1 + 1 \cdot 0.3 + 4 \cdot 0.4 + 9 \cdot 0.2\\ &= 3.7. \end{split} \end{split}$

So, we have:

$\mathbb{V}[X] = \mathbb{E}[X^2] - \mu^2 = 3.7 - 1.7^2 = 0.81.$

Here is how we can find it with Python:

import numpy as np
# The values X can take
xs = np.arange(4)
print('X values: ', xs)
# The probability for each value
ps = np.array([0.1, 0.3, 0.4, 0.2])
print('X probabilities: ', ps)
# And the expectation in a single line
E_X = np.sum(xs * ps)
# The expectation of the square
E_X2 = np.sum(xs ** 2 * ps)
# The variance
V_X = E_X2 - E_X ** 2
print('V[X] = {0:1.2f}'.format(V_X))

X values:  [0 1 2 3]
X probabilities:  [0.1 0.3 0.4 0.2]
V[X] = 0.81


Alternatively, we could use scipy.stats:

X = st.rv_discrete(name='X', values=(xs, ps))
print('V[X] = {0:1.2f}'.format(X.var()))

V[X] = 0.81


The standard deviation is:

print('std of X = {0:1.2f}'.format(X.std()))

std of X = 0.90


Let’s now make a plot. I am going to plot the the PMF of $$X$$ and I am going to mark the position of the expected value along with:

• the expected value minus two standard deviations,

• the expected value plus two standard deviations.

Let’s see what we get:

fig, ax = plt.subplots()
ax.vlines(xs, 0, X.pmf(xs), label='PMF of $X$')
mu = X.expect()
std = X.std()
low = mu - 2 * std
up = mu + 2 * std
ax.plot(mu, 0, 'ro', label='$\mu = \mathbf{E}[X]$')
ax.plot(low, 0, 'gx', label='$\mu - 2\sigma$')
ax.plot(up, 0, 'md', label='$\mu + 2\sigma$')
ax.set_xlabel('$x$')
ax.set_ylabel('$p(x)$')
ax.set_title('Categorical$(0.1, 0.3, 0.4, 0.2)$'.format(theta))
plt.legend(loc='upper left');


We see that, in this case, going two standard deviations below the mean and two standard deviations above the mean captures pretty much all the values.

## Example: Variance of a Binomial random variable¶

Take a Binomial random variable:

$X\sim \text{Binomial}(n, \theta).$

It is not very easy to find the variance of this one. But it is given by the following formula:

$\mathbf{V}[X] = n\theta(1-\theta).$

If you notice, it is basically $$n$$ times the variance of the Bernoulli. This is not accident. The Binomial is actually the sum of $$n$$ independent Bernoulli’s. But we do not know the mathematics to deal with this yet.

Here is how we can get it with scipy.stats:

n = 5
theta = 0.6
X = st.binom(n, theta)
print('E[X] = {0:1.2f}'.format(X.var()))
print('Compare to n * theta = {0:1.2f}'.format(n * theta * (1 - theta)))

E[X] = 1.20
Compare to n * theta = 1.20


Let’s plot the same things we plotted for the categorical:

fig, ax = plt.subplots()
xs = np.arange(n+1)
ax.vlines(xs, 0, X.pmf(xs), label='PMF of $X$')
mu = X.expect()
std = X.std()
low = mu - 2 * std
up = mu + 2 * std
ax.plot(mu, 0, 'ro', label='$\mu = \mathbf{E}[X]$')
ax.plot(low, 0, 'gx', label='$\mu - 2\sigma$')
ax.plot(up, 0, 'md', label='$\mu + 2\sigma$')
ax.set_xlabel('$x$')
ax.set_ylabel('$p(x)$')
ax.set_title(r'Binomial$(n={0:d}, \theta={1:1.2f})$'.format(n, theta))
plt.legend(loc='upper left');


### Questions¶

• Rerun the case of the Binomial with $$n=50$$. Does the shape of the PMF you get look familiar?

## Example: Variance of a Poisson random variable¶

$X\sim \operatorname{Poisson}(\lambda).$

Finding this variance is also non-trival. But it is:

$\mathbf{V}[X] = \lambda.$

Wait a second!!! Didn’t we say that the variance has the square units of $$X$$. If you paid attention the expectation of $$X$$ was also $$\lambda$$. How is it even possible? Well, it is because $$X$$ has no units… It’s just numbers counting events.

Let’s also do it in scipy.stats:

lam = 2.0
X = st.poisson(lam)
print('E[X] = {0:1.2f}'.format(X.var()))

E[X] = 2.00


And let’s visualize everything together like before:

fig, ax = plt.subplots()
xs = np.arange(X.ppf(0.9999)) # I will explain this later
ax.vlines(xs, 0, X.pmf(xs), label='PMF of $X$')
mu = X.expect()
std = X.std()
low = mu - 2 * std
up = mu + 2 * std
ax.plot(mu, 0, 'ro', label='$\mu = \mathbf{E}[X]$')
ax.plot(low, 0, 'gx', label='$\mu - 2\sigma$')
ax.plot(up, 0, 'md', label='$\mu + 2\sigma$')
ax.set_xlabel('$x$')
ax.set_ylabel('$p(x)$')
ax.set_title(r'Poisson$(\lambda={0:1.2f})$'.format(lam))
plt.legend(loc='upper right');


### Question¶

• Rerun the case for the Poisson with a rate parameter $$\lambda = 50$$. Does the shape look familiar?