Probability Distributions Cheat Sheet

from scipy import stats# Normal (Gaussian) distributionnorm = stats.norm(loc=0, scale=1)          # mean=0, std=1print(norm.pdf(0))                          # probability density at x=0print(norm.cdf(1.96))                       # P(X <= 1.96)samples = norm.rvs(size=1000)               # random samples# Binomial distributionbinom = stats.binom(n=10, p=0.5)            # 10 trials, p=0.5 successprint(binom.pmf(5))                         # P(exactly 5 successes)# Poisson distributionpoisson = stats.poisson(mu=3)               # average rate = 3 eventsprint(poisson.pmf(2))                       # P(exactly 2 events)# Uniform distributionuniform = stats.uniform(loc=0, scale=10)    # range [0, 10]

from scipy import statsdata = [23, 25, 21, 22, 24, 20, 26, 27, 19, 24]# Fit a normal distribution to data (maximum likelihood estimate)mu, sigma = stats.norm.fit(data)# Test if data plausibly comes from a normal distributionstat, p_value = stats.shapiro(data)   # Shapiro-Wilk normality testprint(f"p = {p_value:.4f}")           # p < 0.05 suggests non-normal# Kolmogorov-Smirnov test against a specific fitted distributionks_stat, p_value = stats.kstest(data, "norm", args=(mu, sigma))

import numpy as npfrom scipy import statsmean = [0, 0]cov = [[1.0, 0.6], [0.6, 1.0]]  # positive correlation between X and Ymvn = stats.multivariate_normal(mean=mean, cov=cov)sample = mvn.rvs(size=1000, random_state=0)print(mvn.pdf([0.5, 0.5]))          # joint density at a point# Conditional distribution of Y given X=1 (still Gaussian for the MVN family)x_val = 1.0cond_mean = mean[1] + cov[1][0] / cov[0][0] * (x_val - mean[0])cond_var = cov[1][1] - cov[1][0] ** 2 / cov[0][0]cond_dist = stats.norm(loc=cond_mean, scale=np.sqrt(cond_var))# Copula-style dependence: transform correlated normals to arbitrary marginalsu = stats.norm.cdf(sample[:, 0])         # uniform marginal via probability integral transformexp_marginal = stats.expon.ppf(u, scale=2)  # now Exponential, correlation structure preserved

import numpy as npfrom scipy import statsdata = np.random.default_rng(0).gamma(shape=2.0, scale=1.5, size=500)candidates = ["norm", "lognorm", "gamma", "weibull_min", "expon"]results = []for name in candidates:    dist = getattr(stats, name)    params = dist.fit(data)    log_lik = np.sum(dist.logpdf(data, *params))    k = len(params)    aic = 2 * k - 2 * log_lik    ks_stat, ks_p = stats.kstest(data, name, args=params)    results.append((name, aic, ks_p))results.sort(key=lambda r: r[1])   # lower AIC = better fitfor name, aic, ks_p in results:    print(f"{name:12s} AIC={aic:8.1f}  KS p={ks_p:.4f}")

import numpy as npfrom scipy import stats# Transformation of a random variable: if X ~ N(0,1), then Y = X^2 ~ chi-square(1)x = stats.norm.rvs(size=100000, random_state=0)y = x ** 2print(np.mean(y), np.var(y))          # matches chi2(df=1): mean=1, var=2# Log-normal is the exponential transform of a normallog_returns = stats.norm.rvs(loc=0.0005, scale=0.02, size=100000, random_state=1)prices = np.exp(log_returns)          # lognormal-distributed# Order statistics: distribution of the sample max of n uniforms is Beta(n, 1)n = 10maxima = stats.uniform.rvs(size=(100000, n)).max(axis=1)theoretical = stats.beta(a=n, b=1)print(f"empirical mean max: {maxima.mean():.4f}, theoretical: {theoretical.mean():.4f}")

import numpy as npfrom sklearn.mixture import GaussianMixture# Simulate bimodal data: two overlapping clustersrng = np.random.default_rng(0)data = np.concatenate([    rng.normal(-2, 0.8, 300),    rng.normal(3, 1.2, 200),]).reshape(-1, 1)gmm = GaussianMixture(n_components=2, random_state=0).fit(data)print("component means:", gmm.means_.ravel())print("component weights:", gmm.weights_)print("component variances:", gmm.covariances_.ravel())# Use BIC to select the number of componentsbics = [GaussianMixture(n_components=k, random_state=0).fit(data).bic(data) for k in range(1, 6)]best_k = np.argmin(bics) + 1