# Distributions#

This notebook presents various well-defined numeric probability distributions that frequently arise in statistical analysis.

I have written it to make it easy for you to see different distributions and different parameters. Please play with different values to see what happens!

This notebook also demonstrates how to do some plotting with raw matplotlib.

## Math Refresh#

A discrete probability distribution is a distribution over discrete values. These may categorical, ordinal, or numeric. Such a distribution is defined by the probability mass function $$P(x)$$ (for a discrete random variable, $$P(X = x)$$.

A continuous probability distribution is similarly a distribution over continuous values. Because $$P(X = x) = 0$$ for a continuous random variable $$X$$, we instead assign probability mass to intervals. This is done through the distribution function (or cumulative distribution function, or CDF):

$F(x) = P(X \le x)$

Such distributions are often defined in terms of a probability density function $$p(x)$$ such that:

\begin{split}\begin{align*} F(x) & = \int_{-\infty}^x p(x) dx \\ p(x) & = \frac{d}{dx} F(x) \end{align*}\end{split}

The expected value of a random variable is:

$\mathrm{E}[X] = \int x p(x) dx$

For a discrete variable, replace the integral with a sum and the density with the mass:

$\mathrm{E}[X] = \sum_x x P(x)$

Mathematically, we denote that a variable follows a distribution with $$\sim$$. To write that $$X$$ is normally distributed with mean 0 and standard deviation 1, we write:

$X \sim \mathrm{Normal}(0, 1)$

Probability functions often take parameters that govern their behavior. The number, range, and interpretation of parameters varies greatly from function to function.

## Setup#

import numpy as np
import scipy.stats as sps
import seaborn as sns
import matplotlib.pyplot as plt

### Functions#

This notebook is going to be much easier to use (and write) if we define some helper functions.

SciPy provides classes that implement many distributions in a compatible interface. We’re going to use those.

#### Density#

Our first function is going to take a distribution and plot its density. It takes an optional min and max for the $$x$$ axis, and will use a helper function (defined in a moment) to compute default values if they aren’t specified. But the plot function looks like this:

def plot_density(dist, xmin=None, xmax=None, npts=1000, mass=0.9999):
xmin, xmax = find_range(dist, xmin, xmax, mass=mass)
# set up x values
xs = np.linspace(xmin, xmax, npts)
# find corresponding densities
ys = dist.pdf(xs)
# let's start by plotting the mean & median - we plot this first so the density is on top
plt.axvline(dist.median(), linestyle=':', alpha=0.5, label='median')
plt.axvline(dist.mean(), linestyle='--', alpha=0.5, label='mean')
# and plot the density itself
plt.plot(xs, ys)
# and a legend & labels
plt.legend()
plt.xlabel('x')
plt.ylabel('Probability Density')

And plot multiple densities at the same time:

def plot_densities(*dists, xmin=None, xmax=None, npts=1000, mass=0.9999):
ranges = [find_range(d, xmin, xmax, mass) for (l, d) in dists]
xmin = min(r[0] for r in ranges)
xmax = max(r[1] for r in ranges)
# set up x values
xs = np.linspace(xmin, xmax, npts)
# plot y values
for lbl, dist in dists:
# find corresponding densities
ys = dist.pdf(xs)
plt.plot(xs, ys, label=lbl)

# and a legend & labels
plt.legend()
plt.xlabel('x')
plt.ylabel('Probability Density')

#### Distribution#

We also want to be able to plot CDFs:

def plot_cdf(dist, xmin=None, xmax=None, npts=1000, mass=0.9999):
xmin, xmax = find_range(dist, xmin, xmax, mass=mass)
# set up x values
xs = np.linspace(xmin, xmax, npts)
# find corresponding densities
ys = dist.cdf(xs)
# let's start by plotting the mean & median - we plot this first so the density is on top
plt.axvline(dist.median(), linestyle=':', alpha=0.5, label='median')
plt.axvline(dist.mean(), linestyle='--', alpha=0.5, label='mean')
# and plot the density itself
plt.plot(xs, ys)
# and a legend & labels
plt.legend()
plt.xlabel('x')
plt.ylabel('Cumulative Probability')

#### Density and CDF#

Let’s also build a way to plot both the density and the CDF at the same time.

def plot_pdf_cdf(dist, xmin=None, xmax=None, npts=1000, mass=0.9999):
xmin, xmax = find_range(dist, xmin, xmax, mass=mass)
# set up x values
xs = np.linspace(xmin, xmax, npts)

# set up 2 plots
fig, axs = plt.subplots(2, sharex=True, figsize=(6,8))
pdf, cdf = axs  # extract two axes

# let's start by plotting the mean & median - we plot this first so the density is on top
pdf.axvline(dist.median(), linestyle=':', alpha=0.5, label='median')
cdf.axvline(dist.median(), linestyle=':', alpha=0.5, label='median')
pdf.axvline(dist.mean(), linestyle='--', alpha=0.5, label='mean')
cdf.axvline(dist.mean(), linestyle='--', alpha=0.5, label='mean')
# and plot the density itself
pdf.plot(xs, dist.pdf(xs))
# and the CDF
cdf.plot(xs, dist.cdf(xs))
# and a legend & labels
pdf.legend()
pdf.set_ylabel('Probability Density')
cdf.set_xlabel('x')
cdf.set_ylabel('Cumulative Probability')

#### Continuous Range#

Now let’s define that find_range function. If we already have xmin and/or xmax, they are left as-is, so we can manually specify the range. But if they are not, we will use the inverse CDF (which SciPy calls the percent point function ppf) to find values such that the specified probability mass will be between xmin and xmax.

Note: in this function, we reassign the parameter variables. Python is call-by-value, so this does not change the underlying values; it simply reassigns our local variables.

Note 2: if we passed a list as a parameter, and modified the list, it would change it for the caller. This is not a violation of call-by-value; rather, it reflects that object references are passed by value. If we modify a list, it is changed because there’s just one list object shared between caller and callee; if we assign the variable to contain a reference to a new list, the caller doesn’t get to see that, because while it has the list object, it doesn’t know anything about the variable in which the function stores a refernce to it. Just something to be careful of.

Anyway, here’s the function:

def find_range(dist, xmin, xmax, mass=0.999):
half_excluded_mass = (1 - mass) / 2
if xmin is None:
xmin = dist.ppf(half_excluded_mass)
if xmax is None:
xmax = dist.ppf(1 - half_excluded_mass)
return xmin, xmax

#### Mass#

Finally, let’s plot probability mass for discrete distributions:

def plot_mass(dist, xmin=0, xmax=None, mass=0.9999):
if xmax is None:
xmax = dist.ppf(mass + (1 - mass) * 0.5)
# set up x values
xs = np.arange(0, xmax + 1)
# find corresponding densities
ys = dist.pmf(xs)
# let's start by plotting the mean & median - we plot this first so the density is on top
plt.axvline(dist.median(), linestyle=':', alpha=0.5, label='median')
plt.axvline(dist.mean(), linestyle='--', alpha=0.5, label='mean')
# and plot the density itself
plt.bar(xs, ys, width=0.5)
# and a legend & labels
plt.legend()
plt.xlabel('x')
plt.ylabel('Probability')

## Bernoulli Distribution#

We’ll start with a very simple discrete probability distribution: the Bernoulli distribution.

The binomial distribution is a distribution over two outcomes, such as a coin flip. Conventionally, one outcome is called ‘success’ and the other ‘failure’, because it is used to describe the outcome of a Bernoulli trial, which is an action or even tthat can either succeed or fail.

This distribution has a single parameter, $$\theta$$, that is the probability of success:

\begin{split}\begin{align*} P(\mathrm{S}) & = \theta \\ P(\mathrm{F}) & = 1 - \theta \end{align*}\end{split}

We often encode success and failure numerically as 1 and 0, respectively.

plot_mass(sps.bernoulli(0.5), 0, 1)
plt.xticks([0,1], ['F', 'S'])  # drop excess tick labels & label success & failure
plt.show()

## Binomial Distribution#

The binomial distribution is a discrete distribution that describes the number of successes ($$y$$) in a fixed number of Bernoulli trials. It has two parameters: $$\theta \in [0,1]$$ is the success probability, and $$n$$ is the number of trials. Its probability mass function is:

$P(Y = y | \theta, n) = \binom{n}{y} \theta^y (1-\theta)^{(n-y)}$

The probability mass is distributed like this:

plot_mass(sps.binom(10, 0.5), 0, 10)
plt.show()

One thing we can immediately note: a fair coin ($$\theta = 0.5$$), flipped 10 times, will have exactly 5 heads less than 25% of the time!

If we change the success probability, it shifts left or right:

plot_mass(sps.binom(10, 0.3), 0, 10)
plt.show()

Let’s show a more extreme example:

plot_mass(sps.binom(100, 0.1), 0, 100)
plt.show()

## Poisson distribution#

The Poisson distribution frequently shows up when modeling counts. The Poisson has a single parameter, the rate $$\lambda > 0$$, that is the mean (expected value) of the distribution. In terms of a random process, if there is an event that occurs on average $$\lambda$$ times in a fixed time window, then the number of observed occurrances in such a time window is Poisson-distributed.

It is defined by:

$P(X = n | \lambda) = \frac{\lambda^n e^{-\lambda}}{n!}$

Let’s look at one:

plot_mass(sps.poisson(10))

## Negative Binomial#

The negative binomial distribution is another distribution that often shows up in counts. It plays a similar role as the Poisson, but has two parameters so that two negative binomials can have the same mean but different variances.

The parameters are $$r > 0$$ and $$p \in [0,1]$$, and it looks like this:

plot_mass(sps.nbinom(10, 0.5))

In terms of a random process, it’s the number of failures you will see if you run a Bernoilli trial with probability $$p$$ until you see $$r$$ successes. In the distribution shown above, with $$p=0.5$$, it’s the number of of tails you will see if you flip a coin until you get 10 heads.

## Uniform Distribution#

The uniform distribution distributes its probability mass equally across a range. For a discrete variable, every value has equal mass; for a continuous variable, the entire range has equal density.

In order to be defined, it requires both upper and lower limits to the variable.

Here’s a continuous uniform density over $$[0,10]$$:

plot_density(sps.uniform(0, 10))

And the CDF:

plot_cdf(sps.uniform(0, 10))

Straight line! About as simple as it gets.

## Normal Distribution#

A very common continuous distribution is the normal distribution, also known as the Gaussian distribution and the bell curve. It has two parameters, the mean $$\mu$$ and standard deviation $$\sigma$$.

The density of the standard normal ($$\mu=0$$, $$\sigma=1$$) looks like this:

std_norm = sps.norm()
plot_density(std_norm, mass=0.9999)

And its distribution function:

plot_cdf(std_norm, mass=0.9999)

The normal distribution is a location-scale distribution. Such a distribution has location and scale parameters ($$\mu$$ is the location and $$\sigma$$ is the scale), and the distribution can be adjusted through adding the location and multiplying by the scale.

If $$X \sim \mathrm{Normal}(0, 1)$$, and $$Y \sim \mathrm{Normal}(-3, 5)$$, then we can simulate $$y$$ by drawing $$x$$ and computing $$y = 5 - 3 x$$.

Let’s see the density of $$Y$$:

plot_density(sps.norm(-3, 5), mass=0.9999)

The shape is exactly the same; we have just shifted the mean (to -3) and expanded the range of values (from a standard deviation of 1 to 5). This is what happens with any scale-location distribution!

The normal distribution is unbounded ($$p_{\mathrm{Normal}}$$ is defined over the entire real line).

Let’s show off our combined plot, so I can put it in the slide:

plot_pdf_cdf(sps.norm())

## Exponential Distribution#

Another common continuous distribution is the exponential. It is a strongly right-skewed distribution over non-negative numbers defined by a rate parameter $$\lambda > 0$$. Its density is:

$p(x | \lambda) = \lambda e^{-\lambda x}$

Here’s what its density looks like:

plot_density(sps.expon(scale=1))

The basic shape of of the exponential distribution does not change as you change the rate - it just rescales the range over which the curve is distributed, and the $$y$$ values

plot_density(sps.expon(scale=0.5))

Let’s compare them, directly, though:

plot_densities(('λ=1.0', sps.expon(scale=1.0)), ('λ=0.1', sps.expon(scale=0.5)))

## Wrapup#

There are many probability distributions that cover a variety of different processes and models. I have only touched on a few of them here.