Probability and Distributions in Julia

Overview

This tutorial is a reference for the probability you need to do Monte Carlo analysis in this course. It is written for students who have taken a probability or statistics course but may not have used it in a while, and it assumes no prior experience with Distributions.jl.

Everything here is covered or used in lecture, but lecture moves quickly. Come back to this when you need to look something up.

If you want to plot distributions, see the “Plotting Distributions” section of Julia Plotting. This tutorial is about computing with them.

Setup

One package does almost all of the work.

using Distributions
using Random
using Statistics

Random Variables and Distributions

A random variable is a quantity whose value we do not know. A distribution describes how plausible its possible values are.

The distinction that matters in practice is whether the variable is discrete or continuous.

Discrete variables and probability mass

A discrete variable takes values you can list: the number of storms in a year, the number of treatment plants that fail. Its distribution is described by a probability mass function (PMF), which gives the probability of each individual value:

\[p(x) = \mathbb{P}(X = x)\]

These probabilities are genuinely probabilities — they are between 0 and 1, and they sum to 1.

storms = Poisson(3.0)      # average of 3 storms per year
pdf(storms, 2)             # probability of exactly 2 storms
0.22404180765538775

Note that Distributions.jl uses pdf for both cases. For a discrete distribution it returns the mass; for a continuous one it returns the density.

Continuous variables and probability density

A continuous variable takes any value in a range: a concentration, a flow rate, a temperature. Here the probability of any exact value is zero, so we work with a probability density function (PDF) instead:

\[\mathbb{P}(a < X < b) = \int_a^b f(x)\,dx\]

A density is not a probability. It can be larger than 1. Only the area under it is a probability.

load = Normal(9.0, 1.5)    # CBOD load, mg/L
pdf(load, 9.0)             # density at the mean -- not a probability
0.2659615202676218

Narrow the distribution and the density climbs straight past 1, which no probability can do:

pdf(Normal(9.0, 0.2), 9.0)
1.9947114020071635

The Cumulative Distribution Function

The CDF answers the question you usually actually have: how much probability is at or below some value?

\[F(x) = \mathbb{P}(X \le x)\]

cdf(load, 11.0)            # probability the load is at or below 11 mg/L
0.9087887802741321

It works the same way for discrete distributions:

cdf(storms, 2)             # probability of 2 or fewer storms
0.42319008112684353

Exceedance probabilities

Standards are usually written in terms of being worse than a threshold, which is the complement of the CDF:

\[\mathbb{P}(X > x) = 1 - F(x)\]

You can write this directly, but ccdf (“complementary CDF”) is more accurate for small probabilities, because it avoids subtracting two nearly-equal numbers:

1 - cdf(load, 13.0)        # works
0.003830380567589775
ccdf(load, 13.0)           # better
0.003830380567589736

Quantiles

The quantile function is the CDF run backwards. Give it a probability, it returns the value that sits there:

quantile(load, 0.5)        # the median
9.0
quantile(load, 0.95)       # the value exceeded 5% of the time
11.467280440427208

This is how you get the endpoints of an interval. For a 95% interval, take the 2.5% and 97.5% quantiles:

(quantile(load, 0.025), quantile(load, 0.975))
(6.060054023189911, 11.939945976810087)

Summarizing a Distribution

The two summaries you will use constantly:

mean(load), var(load), std(load)
(9.0, 2.25, 1.5)

The average of a function is not the function of the average

This one causes more trouble than anything else on this page.

If you push a random variable through a model \(h\), the average output is generally not the output at the average input:

\[\mathbb{E}[h(X)] \neq h(\mathbb{E}[X])\]

They are equal only when \(h\) is linear. Here is a quick demonstration with \(h(x) = x^2\):

Random.seed!(4750)
draws = rand(load, 100_000)

mean(draws .^ 2), mean(draws)^2
(83.25567014845863, 81.0165967503216)

The two differ by the variance, and for a model as nonlinear as a dissolved oxygen sag curve the gap can be much larger. This is the reason we sample instead of running a model once at the mean.

A Short Catalogue

You do not need to memorize these. You do need to be able to say why you picked one.

Distribution Use it when Constructor
Normal A quantity varies symmetrically about a typical value Normal(μ, σ)
LogNormal A positive quantity varies multiplicatively; right-skewed LogNormal(μ, σ)
Uniform Every value in a range is equally plausible Uniform(a, b)
Exponential Waiting time between independent events Exponential(θ)
Poisson Count of independent events in a fixed window Poisson(λ)
Binomial Number of successes in a fixed number of trials Binomial(n, p)
Beta A proportion, bounded between 0 and 1 Beta(α, β)
Gamma A positive, right-skewed quantity such as a sum of waiting times Gamma(α, θ)

Two cautions worth internalizing.

LogNormal takes the parameters of the underlying normal, not the mean and standard deviation of the distribution itself. LogNormal(log(2), 0.5) has a median of 2, not a mean of 2.

A Normal is unbounded. If you use it for a concentration, it will eventually hand you a negative one. Either check, or use a distribution that cannot:

Random.seed!(4750)
minimum(rand(Normal(2.0, 1.0), 10_000))    # negative concentrations
-1.5881740321384474
truncated(Normal(2.0, 1.0); lower=0.0)     # one way to fix it
Truncated(Distributions.Normal{Float64}(μ=2.0, σ=1.0); lower=0.0)

Drawing Samples

rand draws from any distribution. One value, or many:

Random.seed!(4750)
rand(load)
9.741525584287155
Random.seed!(4750)
rand(load, 5)
5-element Vector{Float64}:
 9.741525584287155
 8.551739787702976
 9.881899367462028
 6.926693118258402
 9.293861462481525

This is the operation at the heart of Monte Carlo: draw inputs, push each one through the model, and look at the distribution of outputs.

Random.seed!(4750)
samples = rand(load, 10_000)
mean(samples), std(samples)
(9.008633740674425, 1.4944133097855608)

Compare those to the true mean(load) and std(load) above. They are close but not exact, and that gap is the Monte Carlo error.

Estimating a probability

To estimate the probability of an event, count the fraction of samples for which it happens:

mean(samples .> 11.0)        # estimated P(load > 11)
0.0886
ccdf(load, 11.0)             # the exact answer, for comparison
0.09121121972586788

Taking the mean of a vector of true/false values gives the fraction that are true. That idiom appears throughout the course.

Seeds and reproducibility

rand uses a pseudorandom number generator: the sequence looks random but is completely determined by a starting seed. Setting the seed makes your results reproducible, by you and by whoever is grading them.

Random.seed!(1)
a = rand(load, 3)
Random.seed!(1)
b = rand(load, 3)
a == b
true

Set a seed in every assignment that uses randomness. It is worth a point on the Monte Carlo rubric, and without it nobody — including you — can reproduce what you did.

Conditional Probability

The probability of \(A\) given that \(B\) happened:

\[\mathbb{P}(A \mid B) = \frac{\mathbb{P}(A \cap B)}{\mathbb{P}(B)}\]

Two events are independent when conditioning changes nothing, so \(\mathbb{P}(A \mid B) = \mathbb{P}(A)\) and \(\mathbb{P}(A \cap B) = \mathbb{P}(A)\mathbb{P}(B)\).

This matters in this course mostly through what it lets you multiply. If annual maxima are independent from year to year, the probability of no exceedance in \(n\) years is \((1-p)^n\), so the probability of at least one is

\[1 - (1 - p)^n.\]

p = 0.01                     # annual exceedance probability, a 100-year event
n = 30                       # design life, years
1 - (1 - p)^n
0.2602996266117198

A “100-year event” has better than a one-in-four chance of showing up during a 30-year design life. Independence is doing real work in that calculation — if exceedances cluster, the answer changes.

Getting Help