Specialized Distributions

Chapter 1: Softmax / Gibbs Distribution

You have a neural network that outputs raw scores (logits) for 5 classes: [2.0, 1.0, 0.1, -1.0, 3.5]. Which class should it pick? The scores aren't probabilities — they can be negative, they don't sum to 1. You need a way to convert arbitrary real numbers into a valid probability distribution.

The softmax function (also called the Gibbs distribution in statistical mechanics) does exactly this. For a vector of scores z = [z₁, ..., z_K], it produces:

The parameter τ (tau) is the temperature. At high temperature (τ → ∞), the distribution becomes uniform — all classes equally likely. At low temperature (τ → 0), all probability concentrates on the highest-scoring class. Think of τ as a "sharpness knob."

Properties

Worked Example

Logits: z = [2.0, 1.0, 0.1]. With τ = 1:

exp(2.0) = 7.389, exp(1.0) = 2.718, exp(0.1) = 1.105. Sum = 11.212.

p = [0.659, 0.242, 0.099]. The highest logit gets 66% of the probability mass.

python
import numpy as np

def softmax(z, tau=1.0):
    e = np.exp(z / tau - np.max(z / tau))  # subtract max for numerical stability
    return e / e.sum()

logits = np.array([2.0, 1.0, 0.1, -1.0, 3.5])
print(softmax(logits, tau=1.0))   # peaked at class 4
print(softmax(logits, tau=0.1))   # nearly one-hot
print(softmax(logits, tau=10.0))  # nearly uniform

Chapter 2: SO(3)/SE(3) Distributions

A robot arm grasps a cup. To describe the cup's orientation, you need three angles (roll, pitch, yaw). But these angles don't live on a flat Euclidean space — they live on a curved surface called a manifold. Specifically, 3D rotations form the group SO(3) (Special Orthogonal group in 3 dimensions).

Why can't we just use a Gaussian on Euler angles? Because of gimbal lock (singularities where two axes align) and wrapping (an angle of 361° is the same as 1°). A Gaussian centered at 350° with σ = 20° would naively split probability between 330-370, but 370° = 10°, so the distribution should wrap around.

The standard solution uses quaternions (4D unit vectors) or rotation matrices (3×3 orthogonal matrices). On these representations, we define distributions directly on the manifold.

Key Distributions on SO(3)

The Bingham distribution is the analog of a Gaussian on the unit sphere. For a unit quaternion q:

where M is a symmetric matrix controlling the concentration and preferred direction. The Matrix Fisher distribution works directly on rotation matrices R:

where F encodes the mean rotation and concentration.

python
import numpy as np
from scipy.spatial.transform import Rotation

def sample_rotation_vmf(mean_axis, kappa, n=100):
    """Sample rotations concentrated around an axis using von Mises-Fisher."""
    # Sample angle from vMF on the circle
    angles = np.random.vonmises(0, kappa, n)
    # Create rotation objects
    rotations = [Rotation.from_rotvec(a * mean_axis) for a in angles]
    return rotations

samples = sample_rotation_vmf(np.array([0,0,1]), kappa=10)

Chapter 3: Wrapped Cauchy

You're studying wind directions at a weather station. Wind can blow from any compass heading: 0° to 360°, and then it wraps back to 0°. You need a probability distribution on the circle. The von Mises distribution (the "circular Gaussian") is one option, but what if your data has heavy tails — occasional extreme deviations from the mean direction?

The Wrapped Cauchy distribution is the circular analog of the Cauchy distribution. Take a Cauchy distribution on the real line, wrap it around the unit circle, and you get a distribution with much heavier tails than von Mises. Its PDF on the circle [0, 2π) is:

Here μ is the mean direction and ρ ∈ [0, 1) is the concentration parameter. When ρ = 0, the distribution is uniform on the circle. As ρ → 1, it concentrates around μ — but always with heavier tails than von Mises at the same concentration.

python
import numpy as np

def wrapped_cauchy_pdf(theta, mu, rho):
    """Wrapped Cauchy PDF on [0, 2pi)."""
    return (1-rho**2) / (2*np.pi * (1+rho**2 - 2*rho*np.cos(theta-mu)))

# Sample via inverse CDF
def sample_wrapped_cauchy(mu, rho, n=1000):
    u = np.random.uniform(0, 1, n)
    return (mu + 2*np.arctan2((1-rho)*np.tan(np.pi*(u-0.5)), 1+rho)) % (2*np.pi)

Chapter 4: Generalized Pareto

You're monitoring a river's water level. Most days it's between 2-4 meters. But during a flood, it could hit 10, 15, even 20 meters. Insurance companies, engineers, and risk analysts don't care about typical days — they need to model the extreme events. How high could the river get in a once-per-century flood?

The Generalized Pareto Distribution (GPD) is the theoretically justified distribution for modeling exceedances over a high threshold. By the Pickands-Balkema-de Haan theorem, the distribution of values above a high threshold converges to a GPD, regardless of the original distribution.

The shape parameter ξ (xi) controls the tail behavior. When ξ > 0, the tail is heavy (power-law decay — extreme events are relatively likely). When ξ = 0, it reduces to the exponential distribution. When ξ < 0, the distribution has a finite upper bound at -σ/ξ.

python
from scipy import stats

# Fit GPD to exceedances above a threshold
data = np.random.pareto(2, 10000)  # heavy-tailed data
threshold = np.percentile(data, 95)
exceedances = data[data > threshold] - threshold

shape, loc, scale = stats.genpareto.fit(exceedances, floc=0)
print(f"Shape xi = {shape:.3f}, Scale sigma = {scale:.3f}")

Chapter 5: Normalizing Flows

What if you want a distribution that can be any shape? Not Gaussian, not Pareto, not anything from a textbook — but a complex, multimodal, asymmetric blob learned from data? Normalizing flows let you start with a simple distribution (like a Gaussian) and warp it through a sequence of invertible transformations until it matches your target.

The key idea: if z ~ N(0, I) and x = f(z) where f is an invertible, differentiable function, then the density of x is given by the change of variables formula:

where J is the Jacobian matrix of partial derivatives. The determinant tells you how much the transformation stretches or compresses space locally.

By stacking K invertible layers, x = f_K ˆ ... ˆ f₁(z), you can learn arbitrarily complex distributions. Each layer must be (1) invertible (so you can compute f^-1), and (2) have a tractable Jacobian determinant (so you can evaluate the density).

python
import torch
from torch import nn

class AffineFlow(nn.Module):
    """Simplest normalizing flow: x = z * exp(s) + t"""
    def __init__(self):
        super().__init__()
        self.s = nn.Parameter(torch.zeros(1))
        self.t = nn.Parameter(torch.zeros(1))

    def forward(self, z):
        x = z * torch.exp(self.s) + self.t
        log_det = self.s  # log|det J| = s
        return x, log_det

    def inverse(self, x):
        z = (x - self.t) * torch.exp(-self.s)
        return z

Chapter 6: Squashed Gaussian

You're training a robot to control a motor. The action (torque) must be between -1 and +1. Your policy network outputs a mean μ and standard deviation σ, defining a Gaussian. But a Gaussian has infinite support — it assigns probability to torques of +50 or -100, which are physically impossible.

The squashed Gaussian fixes this by passing the Gaussian sample through tanh, which squishes any real number into (-1, 1). If u ~ N(μ, σ²) and a = tanh(u), then a is guaranteed to be in (-1, 1). But we need the density of a, not u.

Since tanh'(u) = 1 - tanh²(u), we divide by this factor. In log-probability (which is what we use for policy gradients):

python
import torch
from torch.distributions import Normal

def squashed_gaussian_sample(mu, sigma):
    dist = Normal(mu, sigma)
    u = dist.rsample()            # reparameterized sample
    a = torch.tanh(u)             # squash to (-1, 1)
    # Log-prob with correction
    log_prob = dist.log_prob(u) - torch.log(1 - a**2 + 1e-6)
    return a, log_prob.sum(-1)

Chapter 7: Skew-Normal

Stock returns are famously asymmetric. Losses tend to be larger and more sudden than gains. If you fit a Gaussian, you miss this asymmetry entirely — the bell curve is perfectly symmetric. You need a distribution that can lean to one side.

The Skew-Normal distribution adds a single parameter α to the Gaussian that controls asymmetry. Its PDF is:

where φ(x) is the standard normal PDF and Φ(x) is the standard normal CDF. The parameter α controls skewness: α = 0 gives the standard Gaussian, α > 0 skews right (long right tail), α < 0 skews left.

python
from scipy import stats

# skewnorm(a) where a = alpha (shape parameter)
x = np.linspace(-4, 8, 200)
for alpha in [-4, 0, 4, 8]:
    y = stats.skewnorm.pdf(x, alpha)
    print(f"alpha={alpha}: mean={stats.skewnorm.mean(alpha):.3f}, "
          f"skew={stats.skewnorm.stats(alpha, moments='s')[0]:.3f}")

Chapter 8: Generalized Normal

The Gaussian penalizes outliers with squared distance (x - μ)². The Laplace distribution uses absolute distance |x - μ|, which is gentler on outliers. What if you want something in between? Or even sharper than Gaussian? The Generalized Normal Distribution (GND) provides a single "shape dial."

The parameter β controls the tail weight:

python
from scipy import stats

# scipy.stats.gennorm(beta) is the generalized normal
x = np.linspace(-5, 5, 300)
for beta in [0.5, 1, 2, 5]:
    y = stats.gennorm.pdf(x, beta)
    print(f"beta={beta}: kurtosis={stats.gennorm.stats(beta, moments='k')[0]:.2f}")

Chapter 9: Copulas

Two stocks both have heavy-tailed returns (modeled by, say, Student-t distributions). You know each stock's marginal distribution. But how do they move together? In calm markets they're weakly correlated, but during crashes they drop in lockstep. The dependency structure is separate from the individual distributions.

A copula separates the problem of modeling a multivariate distribution into two parts: (1) the marginals (each variable's individual distribution) and (2) the dependency structure (how they relate). Sklar's theorem guarantees this decomposition always exists:

where F_X, F_Y are the marginal CDFs and C is the copula function mapping [0,1]² → [0,1]. The copula C captures all the dependency information.

Common Copulas

python
from scipy.stats import norm
import numpy as np

def gaussian_copula_sample(rho, n=500):
    """Sample from a Gaussian copula with correlation rho."""
    mean = [0, 0]
    cov = [[1, rho], [rho, 1]]
    z = np.random.multivariate_normal(mean, cov, n)
    u = norm.cdf(z)  # transform to uniform marginals
    return u

Chapter 10: Hawkes Process

An earthquake strikes. Over the next few hours, the rate of aftershocks spikes dramatically, then slowly decays. Each aftershock can itself trigger more aftershocks. A tweet goes viral: every retweet increases the chance of more retweets. These are self-exciting processes — events make future events more likely.

The Hawkes process is a point process where the event rate (intensity) at time t depends on past events. Its conditional intensity is:

Three parameters tell the whole story: μ is the background rate (events that happen spontaneously). α is the excitation (how much each event boosts the rate). β is the decay rate (how fast the excitement fades). The process is stable (doesn't explode) when α < 1.

Applications

python
def simulate_hawkes(mu, alpha, beta, T=100):
    """Simulate a univariate Hawkes process via thinning."""
    events = []
    t = 0
    while t < T:
        lam_bar = mu + alpha * beta * sum(
            np.exp(-beta*(t - ti)) for ti in events)
        dt = np.random.exponential(1 / lam_bar)
        t += dt
        if t >= T: break
        lam_t = mu + alpha * beta * sum(
            np.exp(-beta*(t - ti)) for ti in events)
        if np.random.uniform() < lam_t / lam_bar:
            events.append(t)
    return events

Chapter 11: Lévy / Stable Distributions

The Central Limit Theorem says that if you average many independent, finite-variance random variables, you get a Gaussian. But what if the variance is infinite? City sizes, wealth distributions, earthquake magnitudes, and stock market crashes all have such heavy tails that their variance doesn't converge to a finite number. The Gaussian CLT doesn't apply.

The stable distributions (also called Lévy alpha-stable) generalize the Gaussian to handle these cases. They satisfy a generalized CLT: the sum of heavy-tailed i.i.d. variables converges to a stable distribution.

A stable distribution has four parameters:

When α = 2, you get the Gaussian. When α = 1 and β = 0, you get the Cauchy. For α < 2, the distribution has infinite variance. For α ≤ 1, even the mean is undefined.

python
from scipy import stats

# scipy.stats.levy_stable(alpha, beta, loc, scale)
# alpha=2 → Gaussian, alpha=1, beta=0 → Cauchy
x = np.linspace(-10, 10, 500)
for alpha in [0.5, 1.0, 1.5, 2.0]:
    y = stats.levy_stable.pdf(x, alpha, beta=0)
    # Note: for alpha < 2, variance = infinity

Chapter 12: Diffusion / Score-Based

You want to generate photorealistic images. One approach: start with pure noise and gradually denoise it into a crisp image. But to denoise, you need to know the score function — the gradient of the log-density ∇_x log p(x) — at every noise level. This is the core idea behind diffusion models.

The forward process gradually adds Gaussian noise to data x₀ over T steps:

After enough steps, x_T ≈ N(0, I) — pure noise. The noise schedule {β_t} controls how fast the data is destroyed. Common schedules include:

The reverse process learns to denoise: a neural network s_θ(x_t, t) predicts the score (noise direction) at each step. The score matching loss is:

python
import torch

def linear_schedule(T=1000, start=1e-4, end=0.02):
    return torch.linspace(start, end, T)

def forward_diffusion(x0, t, betas):
    alphas = 1 - betas
    alpha_bar = torch.cumprod(alphas, dim=0)[t]
    noise = torch.randn_like(x0)
    xt = torch.sqrt(alpha_bar) * x0 + torch.sqrt(1 - alpha_bar) * noise
    return xt, noise

Chapter 13: Taxonomy & Connections

We've covered 12 specialized distributions. Each exists because some domain has structure that standard distributions can't capture. Here's the complete taxonomy:

By Domain Structure

Relationship Map

When to Reach for Each

Connections to Other Lessons

These distributions appear throughout modern ML and robotics:

"All models are wrong, but some are useful." — George Box