用Python和NumPy手把手教你计算多元正态分布的概率（附身高体重实例）

用Python和NumPy手把手实现多元正态分布概率计算

假设你正在分析一组包含身高体重数据的二维数据集，需要计算特定样本点的概率密度。传统统计学教材中复杂的矩阵运算公式常常让人望而生畏，而今天我们将用Python代码将这些数学符号转化为可执行的程序逻辑。

1. 理解多元正态分布的核心要素

多元正态分布是单变量正态分布在多维空间的自然延伸。想象一个三维空间中的钟形曲面，这个曲面在任意方向上的切片都是一个正态分布曲线。要完整描述这个分布，我们需要两个关键参数：

• 均值向量μ：决定分布中心点的位置
• 协方差矩阵Σ：控制分布的形态和方向

在身高体重数据集中，假设我们观察到：

import numpy as np # 均值向量 [身高均值, 体重均值] mu = np.array([175.0, 71.3]) # 协方差矩阵 [[身高方差, 协方差], [协方差, 体重方差]] sigma = np.array([[874, 327], [327, 929]])

2. 概率密度函数的代码实现

多元正态分布的概率密度函数(PDF)公式看似复杂，但可以分解为几个可计算的组成部分：

def multivariate_normal_pdf(x, mu, sigma): """ 计算多元正态分布的概率密度函数 参数: x: 样本点 (n维向量) mu: 均值向量 sigma: 协方差矩阵 返回: 概率密度值 """ # 维度检查 d = len(mu) if x.shape != (d,) or sigma.shape != (d,d): raise ValueError("维度不匹配") # 计算核心组成部分 diff = x - mu inv_sigma = np.linalg.inv(sigma) det_sigma = np.linalg.det(sigma) # 指数部分计算 exponent = -0.5 * np.dot(np.dot(diff.T, inv_sigma), diff) # 归一化系数 normalization = 1.0 / ((2 * np.pi) ** (d/2) * np.sqrt(det_sigma)) return normalization * np.exp(exponent)

让我们拆解这个实现的关键步骤：

1. 协方差矩阵求逆：np.linalg.inv()计算Σ⁻¹
2. 行列式计算：np.linalg.det()获取|Σ|
3. 马氏距离计算：(x-μ)ᵀΣ⁻¹(x-μ)衡量样本与均值的距离
4. 归一化处理：确保概率密度积分为1

3. 实际应用：身高体重数据分析

假设我们有一个身高178cm，体重73kg的个体，计算其概率密度：

# 样本点 [身高, 体重] sample = np.array([178, 73]) # 计算概率密度 pdf_value = multivariate_normal_pdf(sample, mu, sigma) print(f"概率密度值: {pdf_value:.6f}")

典型输出可能类似于：

概率密度值: 0.000432

注意：概率密度值本身不是概率，只有在积分区间内才有概率意义。对于连续分布，单点的概率理论为0。

4. 可视化与结果验证

为了验证我们的实现是否正确，可以对比SciPy库的内置函数：

from scipy.stats import multivariate_normal # 使用SciPy验证 rv = multivariate_normal(mu, sigma) scipy_pdf = rv.pdf(sample) print(f"我们的实现: {pdf_value:.6f}") print(f"SciPy的结果: {scipy_pdf:.6f}") print(f"差异: {abs(pdf_value - scipy_pdf):.2e}")

输出差异应该在数值计算的误差范围内（通常<1e-10）。

5. 性能优化技巧

当处理高维数据或大量计算时，可以考虑以下优化：

1. 矩阵分解预计算：

# Cholesky分解提高求逆效率 L = np.linalg.cholesky(sigma) inv_sigma = np.linalg.solve(L.T, np.linalg.solve(L, np.eye(2)))

2. 对数概率计算：

def log_multivariate_normal_pdf(x, mu, sigma): diff = x - mu L = np.linalg.cholesky(sigma) log_det = 2 * np.sum(np.log(np.diag(L))) sol = np.linalg.solve(L, diff) return -0.5 * (len(mu)*np.log(2*np.pi) + log_det + np.dot(sol.T, sol))

3. 批量计算：

def batch_multivariate_normal_pdf(X, mu, sigma): """ X: (n_samples, n_dims) """ diff = X - mu inv_sigma = np.linalg.inv(sigma) exponent = -0.5 * np.sum(diff @ inv_sigma * diff, axis=1) normalization = 1.0 / ((2 * np.pi) ** (len(mu)/2) * np.sqrt(np.linalg.det(sigma))) return normalization * np.exp(exponent)

6. 实际案例：异常检测应用

多元正态分布常用于异常检测。假设我们有以下数据：

我们可以计算每个样本的概率密度，并设定阈值：

samples = np.array([[175, 71], [182, 68], [190, 100], [165, 90]]) probs = batch_multivariate_normal_pdf(samples, mu, sigma) threshold = 0.0001 # 通过历史数据确定 for i, (sample, prob) in enumerate(zip(samples, probs)): status = "异常" if prob < threshold else "正常" print(f"样本{i+1}: {sample}, 概率密度: {prob:.2e}, 判断: {status}")

7. 常见问题排查

在实现过程中可能会遇到以下问题：

1. 协方差矩阵不正定：

   • 症状：np.linalg.cholesky抛出LinAlgError
   • 解决：添加小的正则项sigma += 1e-6 * np.eye(len(mu))

2. 概率密度计算为0：

   • 检查是否数值下溢，考虑使用对数概率

3. 维度不匹配错误：

   • 确保所有向量和矩阵的维度一致
   • 使用np.atleast_2d和np.squeeze进行维度调整

4. 性能瓶颈：

   • 对于高维数据，考虑使用对角协方差矩阵
   • 使用np.einsum优化矩阵运算