Lecture 9: randomness

Goal: derive a mathematical representation for “random” variables and their statistics.

References

Section 2.4: Stengel, R. F. Optimal Control and Estimation. (Dover, New York, 1994).

Philosophy

Regardless of whether the universe is “random”, we can use mathematical models of “randomness” in systems to represent unmodeled phenomena (e.g. vibrations, deformations, friction, backlash) or uncertainty (e.g. parameters, estimates).

Definitions

Definition: a random variable is a function $x:\Omega\rightarrow\R$ over a sample space $\Omega$ with associated probabilities.

We will deal with probability distributions over discrete and continuous spaces.

discrete probabilities	continuous probabilities
$\#(\Omega) < \infty$	$\#(\Omega) = \infty$, e.g. $\Omega = \R^n$
each $\omega\in\Omega$ has probability $P(\omega)\in[0,1]$ and $\sum_{\omega\in\Omega} P(\omega) = 1$	each $\omega\in\Omega$ has density $\rho(\omega)\in[0,\infty)$ and $\int_\Omega \rho(\omega) d\omega = 1$.
each event $W\subset\Omega$ has probability $P(W) = \sum_{\omega\in\Omega} P(\omega)$	each event $W \subset\Omega$ has probability $P(W) = \int_W \rho(\omega) d\omega$

It is important to know that rigorously formalizing continuous probabilities requires more work to avoid fundamental problems like the Banach-Tarski paradox. The rigorous formalization is called measure theory. But in this course we will only restrict continuous probability distributions to Gaussians, which avoids these technical issues.

Example (discrete probabilities): consider the random process of flipping a coin. Colloquially, we understand there are two outcomes: $\Omega = \{ H, T \}$. If it’s a fair coin then $P(H) = P(T) = 1/2$. An example random variable could be the “payout” of a gambling game, e.g. $x:\Omega\rightarrow\R$ defined by $x(H) = \$1$ and $x(T) = -\$2$.

Example (continuous probabilities): consider a Gaussian (also termed normal) probability distribution with mean $\mu\in\R^n$ and covariance $\Sigma\in\R^{n\times n}$ denoted $x \sim \mathcal{N}(\mu,\Sigma)$. Then the sample space is $\Omega = \R^n$ and for all $W\subset\Omega$ we have

\[ P(W) = \left[ (2\pi)^n \det\Sigma \right]^{-1/2} \int_W \exp\left( -\frac{1}{2} (\omega - \mu)^\top \Sigma^{-1} (\omega - \mu)\right) d\omega. \]

A random variable $x:\Omega\rightarrow\R^m$ for some $m\in\N$ with this probability distribution could be as simple as $x(\omega) = \omega$ (in which case $m = n$) or $x(\omega) = M\omega$ (in which case $M\in\R^{m\times n}$ with $m\in\N$ arbitrary).

Statistics

A statistic is any function of a random variable $x:\Omega\rightarrow\R^n$. There are many possible choices, including order statistics (e.g. median, interquartile). But we will focus on the mean and (co)variance that define the Gaussian distribution.

Definition: the expectation (or mean) $\text{E}[x\in\R^{n}]$ of a random variable $x:\Omega\rightarrow\R^n$ is defined as follows.

\[ \text{discrete probabilities:}\ \text{E}[x] = \sum_{\omega\in\Omega} x(\omega) P(\omega) \] \[ \text{continuous probabilities:}\ \text{E}[x] = \int_{\Omega} x(\omega) \rho(\omega) d\omega \]

Definition: the (co)variance $\text{Cov}[x]\in\R^{n\times n}$ of a random variable $x:\Omega\rightarrow\R^n$ is defined as follows.

\[ \text{Cov}[x] = \text{E}[(x - \text{E}[x])(x - \text{E}[x])^\top]\in\R^{n\times n} \]

Example (Gaussian random variable): if $x\sim\mathcal{N}(\mu,\Sigma)$ then $\text{E}[x] = \mu$ and $\text{Cov}[x] = \Sigma$. The contours of equal probability density for this distribution are nested ellipses with principal axes aligned with eigenvectors of $\Sigma$. To see this, note that $\Sigma^\top = \Sigma$ from the definition of covariance, so $\Sigma$ is diagonalizable by an orthogonal transformation, i.e. $\Sigma$ has an orthonormal basis of eigenvectors $V\in\R^{n\times n}$ and $V\Sigma V^\top = \text{diag}(\sigma_1^2,\sigma_2^2,\dots,\sigma_n^2)$ is the covariance in this basis. The number $\sigma_i^2$ is the variance of the random variable $\nu_i = v_i^\top x$.

This geometry is illustrated in the following plots.

Code

from imports import *

# Define principal axes (normalized)
v1 = np.array([1, 1]) / np.sqrt(2)  # (1,1) normalized
v2 = np.array([-1, 1]) / np.sqrt(2)  # (-1,1) normalized

# Define variances along principal axes
sigma1_sq = 2  # variance along v1
sigma2_sq = 1  # variance along v2

# Construct covariance matrix: Σ = σ1²·v1·v1^T + σ2²·v2·v2^T
Sigma = sigma1_sq * np.outer(v1, v1) + sigma2_sq * np.outer(v2, v2)

# Mean (centered at origin)
mu = np.array([0.0, 0.0])

# Create grid for evaluation
x1_range = np.linspace(-4, 4, 200)
x2_range = np.linspace(-4, 4, 200)
X1, X2 = np.meshgrid(x1_range, x2_range)

# Stack into column vectors for vectorized computation
omega = np.stack([X1.ravel(), X2.ravel()], axis=0)  # 2 x (200*200)
omega_centered = omega - mu[:, np.newaxis]  # center around mean

# Compute Gaussian density: ρ(ω) = (2π)^(-n/2) |det(Σ)|^(-1/2) exp(-½(ω-μ)^T Σ^(-1) (ω-μ))
Sigma_inv = inv(Sigma)
det_Sigma = det(Sigma)
norm_const = (2 * np.pi) ** (-1) * (det_Sigma) ** (-0.5)

# Compute quadratic form: (ω-μ)^T Σ^(-1) (ω-μ) for all points
quad_form = np.sum(omega_centered * (Sigma_inv @ omega_centered), axis=0)

# Compute density
rho = norm_const * np.exp(-0.5 * quad_form)
rho = rho.reshape(X1.shape)

# Visualization 1: 3D plot with contours
fig = plt.figure(figsize=(12, 5))

# 3D surface plot
ax1 = fig.add_subplot(121, projection='3d')
surf = ax1.plot_surface(X1, X2, rho, cmap='viridis', alpha=0.8, linewidth=0, antialiased=True)

# Overlay contour lines on the surface
contour_levels = np.linspace(rho.min(), rho.max(), 10)
ax1.contour(X1, X2, rho, levels=contour_levels, colors='black', alpha=0.3, linewidths=0.5, offset=rho.min())

ax1.set_xlabel(r'$x_1(\omega)$', fontsize=12, labelpad=8)
ax1.set_ylabel(r'$x_2(\omega)$', fontsize=12, labelpad=8)
ax1.set_zlabel(r'$\rho$', fontsize=12, labelpad=8)
ax1.set_title('3D Gaussian Density with Contours')

# Visualization 2: Contour plot with principal axes
ax2 = fig.add_subplot(122)
contour = ax2.contour(X1, X2, rho, levels=15, cmap='viridis')
ax2.clabel(contour, inline=True, fontsize=8)

# Draw principal axes arrows with lengths proportional to variances
# Arrow length = sqrt(variance) for visualization
arrow_scale1 = np.sqrt(sigma1_sq)  # √2 for v1
arrow_scale2 = np.sqrt(sigma2_sq)  # 1 for v2

# Draw v1 arrow
ax2.arrow(0, 0, arrow_scale1 * v1[0], arrow_scale1 * v1[1],
          head_width=0.15, head_length=0.15, fc='red', ec='red', linewidth=2, zorder=5)
ax2.text(arrow_scale1 * v1[0] + 0.2, arrow_scale1 * v1[1] + 0.2, r'$v_1$',
         fontsize=12, color='red', weight='bold')

# Draw v2 arrow
ax2.arrow(0, 0, arrow_scale2 * v2[0], arrow_scale2 * v2[1],
          head_width=0.15, head_length=0.15, fc='blue', ec='blue', linewidth=2, zorder=5)
ax2.text(arrow_scale2 * v2[0] - 0.3, arrow_scale2 * v2[1] + 0.2, r'$v_2$',
         fontsize=12, color='blue', weight='bold')

ax2.set_xlabel(r'$x_1(\omega)$', fontsize=12)
ax2.set_ylabel(r'$x_2(\omega)$', fontsize=12)
ax2.set_title('Contour Plot with Principal Axes')
ax2.set_aspect('equal')
ax2.grid(True, alpha=0.3)

plt.tight_layout()

# Display with alt text
show_with_alt(fig,
    "Two visualizations of a 2D Gaussian probability density. Left: 3D surface plot showing the density "
    "as a height map with contour lines overlaid on the surface. The x and y axes are labeled x₁(ω) and x₂(ω), "
    "and the z-axis is labeled ρ. Right: 2D contour plot of the same density with principal axes shown as arrows. "
    "The red arrow labeled v₁ points in the direction (1,1) with length proportional to variance 2, and the blue "
    "arrow labeled v₂ points in the direction (-1,1) with length proportional to variance 1. Both plots show "
    "the elliptical shape of the Gaussian distribution oriented along the principal axes."
)

Two visualizations of a 2D Gaussian probability density. Left: 3D surface plot showing the density as a height map with contour lines overlaid on the surface. The x and y axes are labeled x₁(ω) and x₂(ω), and the z-axis is labeled ρ. Right: 2D contour plot of the same density with principal axes shown as arrows. The red arrow labeled v₁ points in the direction (1,1) with length proportional to variance 2, and the blue arrow labeled v₂ points in the direction (-1,1) with length proportional to variance 1. Both plots show the elliptical shape of the Gaussian distribution oriented along the principal axes.

Gaussians

A final property of Gaussian random variables that will be critical in subsequent lectures is that an affine transformation of a Gaussian is also Gaussian. This property is very special and does not arise in general.

Fact: if $x\sim\mathcal{N}(\mu,\Sigma)$ and $y = M x + b$ then $y \sim\mathcal{N}(M\mu+b,M\Sigma M^\top)$, i.e. $y$ is a Gaussian random variable and $\text{E}[y] = M\mu + b$, $\text{Cov}[y] = M\Sigma M^\top$.

discrete probabilities	continuous probabilities
\(\#(\Omega) < \infty\)	\(\#(\Omega) = \infty\), e.g. \(\Omega = \R^n\)
each \(\omega\in\Omega\) has probability \(P(\omega)\in[0,1]\) and \(\sum_{\omega\in\Omega} P(\omega) = 1\)	each \(\omega\in\Omega\) has density \(\rho(\omega)\in[0,\infty)\) and \(\int_\Omega \rho(\omega) d\omega = 1\).
each event \(W\subset\Omega\) has probability \(P(W) = \sum_{\omega\in\Omega} P(\omega)\)	each event \(W \subset\Omega\) has probability \(P(W) = \int_W \rho(\omega) d\omega\)