Distributions and Sampling

A continuous random variable $X$ has probability density function $f_X(x)$ satisfying:

$$P(a \le X \le b) = \int_a^b f_X(x)\,dx, \quad f_X(x) \ge 0, \quad \int_{-\infty}^{\infty} f_X(x)\,dx = 1$$

In scipy.stats, every distribution object provides pdf(x), cdf(x), rvs(size), fit(data), and ppf(q) (quantile function).

from scipy.stats import norm
X = norm(loc=0, scale=1)      # standard normal
X.pdf(0.0)                    # f(0) = 0.3989...
X.cdf(1.96)                   # P(X <= 1.96) = 0.975
samples = X.rvs(size=10000)   # draw 10000 samples

The cumulative distribution function (CDF) is:

$$F_X(x) = P(X \le x) = \int_{-\infty}^{x} f_X(t)\,dt$$

Its inverse $F_X^{-1}(q)$ is the quantile function (ppf in SciPy): for $0 < q < 1$, $F_X^{-1}(q)$ returns the value $x$ such that $P(X \le x) = q$.

from scipy.stats import norm
norm.ppf(0.975)    # 1.96 — the 97.5th percentile
norm.isf(0.025)    # 1.96 — inverse survival function

The Gaussian distribution $X \sim \mathcal{N}(\mu, \sigma^2)$ has PDF:

$$f_X(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$

The $Q$-function $Q(x) = 1 - \Phi(x)$ gives the tail probability and appears throughout BER analysis.

from scipy.stats import norm
from scipy.special import erfc
Q = lambda x: 0.5 * erfc(x / np.sqrt(2))
# Equivalent:
Q_scipy = lambda x: norm.sf(x)           # survival function

A circularly symmetric complex Gaussian (CSCG) random variable $Z \sim \mathcal{CN}(0, \sigma^2)$ is defined as:

$$Z = X + jY, \quad X, Y \sim \mathcal{N}(0, \sigma^2/2) \text{ independent}$$

The envelope $|Z|$ follows a Rayleigh distribution and the power $|Z|^2$ follows an exponential distribution with mean $\sigma^2$.

rng = np.random.default_rng(42)
sigma2 = 1.0
z = rng.standard_normal(10000) * np.sqrt(sigma2/2) \
  + 1j * rng.standard_normal(10000) * np.sqrt(sigma2/2)
# |z| is Rayleigh, |z|^2 is exponential

A random vector $\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}, \mathbf{R})$ has joint PDF:

$$f(\mathbf{x}) = \frac{1}{(2\pi)^{n/2} |\mathbf{R}|^{1/2}} \exp\!\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^T \mathbf{R}^{-1} (\mathbf{x}-\boldsymbol{\mu})\right)$$

where $\mathbf{R} = E[(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^T]$ is the covariance matrix (positive semi-definite).

from scipy.stats import multivariate_normal
R = np.array([[1.0, 0.8], [0.8, 1.0]])
rv = multivariate_normal(mean=[0, 0], cov=R)
samples = rv.rvs(size=5000)

To generate $\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}, \mathbf{R})$:

1. Compute the Cholesky factor: $\mathbf{R} = \mathbf{L}\mathbf{L}^T$
2. Draw i.i.d. standard normals: $\mathbf{z} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})$
3. Set $\mathbf{x} = \boldsymbol{\mu} + \mathbf{L}\mathbf{z}$

Then $\mathrm{Cov}[\mathbf{x}] = \mathbf{L} \cdot \mathbf{I} \cdot \mathbf{L}^T = \mathbf{R}$.

L = np.linalg.cholesky(R)    # lower triangular
z = rng.standard_normal((n, N_samples))
x = mu[:, None] + L @ z      # each column is a correlated sample