用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]], size=1000 ) # 估计参数 mu_est = np.mean(train_data, axis=0) 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, projection='3d') ax1.plot_surface(x, y, rv.pdf(pos), cmap='viridis') 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)) # 条件数很大说明矩阵接近奇异
