用Python和NumPy手把手教你计算多元高斯分布的概率密度（附完整代码）

用Python和NumPy手把手实现多元高斯分布概率密度计算第一次看到多元高斯分布的公式时那个包含矩阵运算的复杂表达式确实让人望而生畏。作为机器学习的基础概率模型它描述的是多维空间中数据的分布规律。但别担心今天我们就用Python和NumPy库把这个看似复杂的数学公式转化为可运行的代码。1. 准备工作理解核心概念多元高斯分布Multivariate Gaussian Distribution是单变量高斯分布在多维空间的推广。想象一下我们不再只是描述一个变量的分布比如身高而是同时描述多个相关变量的联合分布比如身高和体重。关键参数解析均值向量μ每个维度的平均值组成的向量协方差矩阵Σ描述各维度之间相关性的对称矩阵举个例子假设我们研究一个国家成年人的身高和体重import numpy as np # 均值向量身高(cm)和体重(kg)的平均值 mu np.array([170, 65]) # 协方差矩阵 sigma np.array([[100, 30], # 身高方差100体重方差25 [30, 25]]) # 协方差30表示正相关2. 概率密度函数拆解实现让我们把多元高斯分布的公式分解为可计算的几个部分f(x) (1/(2π)^(D/2)) * (1/|Σ|^(1/2)) * exp(-1/2 * (x-μ)^T * Σ^(-1) * (x-μ))2.1 计算协方差矩阵的行列式行列式|Σ|反映了数据分布的体积大小。在NumPy中计算很简单def compute_determinant(sigma): return np.linalg.det(sigma) # 示例计算 sigma np.array([[4, 2], [2, 3]]) print(行列式值:, compute_determinant(sigma)) # 输出: 8.02.2 计算协方差矩阵的逆矩阵Σ^(-1)是协方差矩阵的逆用于衡量点与均值之间的马氏距离def compute_inverse(sigma): return np.linalg.inv(sigma) # 示例计算 print(逆矩阵:\n, compute_inverse(sigma))2.3 实现指数部分计算这部分计算点x与均值μ之间的马氏距离def compute_exponent(x, mu, sigma_inv): diff x - mu return -0.5 * np.dot(np.dot(diff.T, sigma_inv), diff) # 示例计算 x np.array([1, 2]) mu np.array([0, 0]) sigma_inv compute_inverse(sigma) print(指数部分:, compute_exponent(x, mu, sigma_inv))3. 完整概率密度函数实现现在我们把所有部分组合起来def multivariate_gaussian(x, mu, sigma): 计算多元高斯分布的概率密度 参数: x: 输入向量(D,) mu: 均值向量(D,) sigma: 协方差矩阵(D,D) 返回: 概率密度值 D len(mu) sigma_inv np.linalg.inv(sigma) det np.linalg.det(sigma) # 常数部分 const 1 / ((2 * np.pi) ** (D/2) * np.sqrt(det)) # 指数部分 diff x - mu exponent -0.5 * np.dot(np.dot(diff.T, sigma_inv), diff) return const * np.exp(exponent)使用示例# 定义参数 mu np.array([0, 0]) sigma np.array([[1, 0.5], [0.5, 1]]) # 计算点[0.5, 0.5]的概率密度 x np.array([0.5, 0.5]) print(概率密度:, multivariate_gaussian(x, mu, sigma))4. 实际应用案例异常检测多元高斯分布在异常检测中非常有用。我们可以计算新数据点的概率密度低于阈值则判定为异常。# 训练数据假设已经标准化 train_data np.random.multivariate_normal( mean[0, 0], cov[[1, 0.5], [0.5, 1]], size1000 ) # 估计参数 mu_est np.mean(train_data, axis0) sigma_est np.cov(train_data.T) # 测试数据 test_point np.array([2.5, 2.5]) # 可能异常 prob multivariate_gaussian(test_point, mu_est, sigma_est) # 设置阈值通常选择训练集概率的某个百分位 threshold 0.01 print(异常 if prob threshold else 正常)5. 性能优化与数值稳定性当维度很高时直接计算可能会遇到数值问题。以下是几个优化技巧5.1 避免直接求逆使用Cholesky分解提高计算效率和稳定性def multivariate_gaussian_optimized(x, mu, sigma): D len(mu) L np.linalg.cholesky(sigma) # Cholesky分解 diff x - mu sol np.linalg.solve(L, diff) exponent -0.5 * np.dot(sol.T, sol) const 1 / ((2 * np.pi) ** (D/2) * np.prod(np.diag(L))) return const * np.exp(exponent)5.2 对数概率计算对于非常小的概率值使用对数空间计算更稳定def log_multivariate_gaussian(x, mu, sigma): D len(mu) sigma_inv np.linalg.inv(sigma) sign, logdet np.linalg.slogdet(sigma) diff x - mu exponent -0.5 * np.dot(np.dot(diff.T, sigma_inv), diff) log_const -0.5 * (D * np.log(2 * np.pi) logdet) return log_const exponent6. 可视化理解可视化能帮助我们直观理解多元高斯分布的形状import matplotlib.pyplot as plt from scipy.stats import multivariate_normal # 创建网格 x, y np.mgrid[-3:3:.1, -3:3:.1] pos np.dstack((x, y)) # 计算概率密度 rv multivariate_normal([0, 0], [[1, 0.5], [0.5, 1]]) fig plt.figure(figsize(10, 5)) ax1 fig.add_subplot(121, projection3d) ax1.plot_surface(x, y, rv.pdf(pos), cmapviridis) ax2 fig.add_subplot(122) ax2.contourf(x, y, rv.pdf(pos)) plt.show()这段代码会生成3D曲面图和等高线图清晰展示概率密度在不同区域的分布情况。7. 常见问题排查在实际实现中你可能会遇到以下问题问题1协方差矩阵不是正定的解决方法添加小的正则项如 sigma 1e-6 * np.eye(D)问题2概率计算结果为0解决方法使用对数概率计算或检查输入数据范围问题3高维情况下计算缓慢解决方法使用Cholesky分解优化或考虑对角协方差矩阵问题4逆矩阵计算不稳定解决方法使用np.linalg.pinv计算伪逆或检查协方差矩阵条件数# 检查协方差矩阵条件数 print(条件数:, np.linalg.cond(sigma)) # 条件数很大说明矩阵接近奇异