用Python和NumPy实战多元高斯分布：从公式推导到身高体重概率计算

在数据科学和机器学习领域，理解如何计算多元高斯分布的概率密度是基础中的基础。想象一下，你手头有一组身高体重的二维数据，已知它们服从某个特定的多元正态分布。现在，你需要计算某个特定身高体重范围（比如身高70-80英寸，体重65-75公斤）出现的概率。这听起来像是一个纯理论问题，但实际上，这正是许多实际应用场景的核心——从异常检测到分类器设计，都离不开这个基础技能。

本文将带你从零开始，用Python和NumPy一步步实现这个计算过程。不同于单纯的理论讲解，我们会聚焦于 如何将数学公式转化为可运行的代码 ，并在最后给出一个完整的、可直接复用的实现方案。无论你是机器学习初学者，还是希望巩固基础的数据科学从业者，这篇文章都将为你提供实用的技术指导。

1. 理解多元高斯分布的核心公式

多元高斯分布的概率密度函数看起来可能有些吓人，但当我们拆解它的各个部分时，会发现每个组件都有明确的物理意义。让我们先看看这个标准公式：

$$ f_{\mu, \Sigma}(x) = \frac{1}{(2\pi)^{D/2}|\Sigma|^{1/2}} \exp\left(-\frac{1}{2}(x-\mu)^T \Sigma^{-1} (x-\mu)\right) $$

这个公式中，每个符号代表什么？它们又如何影响最终的分布形状？

• x ：D维随机变量向量（在我们的身高体重案例中，D=2）
• μ ：均值向量，表示分布的中心位置
• Σ ：协方差矩阵，决定分布的形态和方向
• |Σ| ：协方差矩阵的行列式
• Σ⁻¹ ：协方差矩阵的逆矩阵

为了更好地理解，让我们用Python先定义这些基本组件。假设我们有以下参数（与原始内容中给出的相同）：

import numpy as np

# 定义均值向量和协方差矩阵
mu = np.array([75.0, 71.3])  # 身高和体重的均值
sigma = np.array([[874, 327], [327, 929]])  # 协方差矩阵

1.1 协方差矩阵的几何意义

协方差矩阵Σ决定了多元高斯分布的形状。在我们的二维案例中：

$$ \Sigma = \begin{bmatrix} \sigma_{11} & \sigma_{12} \ \sigma_{21} & \sigma_{22} \end{bmatrix} $$

其中：

• σ₁₁和σ₂₂分别控制身高和体重维度的"胖瘦"（方差）
• σ₁₂=σ₂₁决定两个维度之间的相关性（协方差）

让我们可视化这个协方差矩阵对应的分布形状：

import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

# 创建多元高斯分布实例
rv = multivariate_normal(mu, sigma)

# 生成网格点
x, y = np.mgrid[60:90:0.5, 50:90:0.5]
pos = np.dstack((x, y))

# 绘制等高线图
plt.contourf(x, y, rv.pdf(pos), levels=15, cmap='Purples')
plt.xlabel('Height (inches)')
plt.ylabel('Weight (kg)')
plt.colorbar(label='Probability Density')
plt.title('2D Gaussian Distribution of Height and Weight')
plt.show()

这段代码会生成一个二维高斯分布的等高线图，清晰地展示出身高和体重之间的相关性——当身高增加时，体重也倾向于增加，这正是协方差矩阵中非对角线元素为正所表示的关系。

2. 实现概率密度函数的计算

现在，我们来实现这个概率密度函数。虽然SciPy已经提供了现成的实现，但自己动手写一遍能加深理解。我们将分步骤构建这个函数。

2.1 计算关键组件

首先，我们需要计算几个关键量：

1. 协方差矩阵的行列式（|Σ|）
2. 协方差矩阵的逆（Σ⁻¹）
3. 二次型：(x-μ)ᵀΣ⁻¹(x-μ)

def multivariate_gaussian_pdf(x, mu, sigma):
    """
    计算多元高斯分布的概率密度函数
    参数:
        x: 输入向量(shape: (D,))
        mu: 均值向量(shape: (D,))
        sigma: 协方差矩阵(shape: (D,D))
    返回:
        概率密度值
    """
    D = len(mu)
    x = np.asarray(x)
    
    # 1. 计算行列式
    det = np.linalg.det(sigma)
    
    # 2. 计算逆矩阵
    inv_sigma = np.linalg.inv(sigma)
    
    # 3. 计算归一化常数
    norm_const = 1.0 / ((2 * np.pi) ** (D/2) * np.sqrt(det))
    
    # 4. 计算指数部分
    x_mu = x - mu
    exponent = -0.5 * x_mu.T @ inv_sigma @ x_mu
    
    # 5. 计算最终概率密度
    return norm_const * np.exp(exponent)

让我们测试一下这个函数，计算身高=75英寸，体重=71.3kg处的概率密度（正好是均值点）：

test_point = np.array([75.0, 71.3])
print(f"PDF at mean: {multivariate_gaussian_pdf(test_point, mu, sigma):.6f}")

输出应该与SciPy的实现一致：

from scipy.stats import multivariate_normal
rv = multivariate_normal(mu, sigma)
print(f"Scipy PDF at mean: {rv.pdf(test_point):.6f}")

2.2 处理数值稳定性问题

在实际应用中，我们可能会遇到数值不稳定的情况，特别是当协方差矩阵接近奇异（行列式接近零）时。我们可以通过以下改进增强鲁棒性：

1. 添加一个小常数到对角线元素（正则化）
2. 使用伪逆代替常规逆
3. 使用对数空间计算避免数值下溢

改进后的版本：

def safe_multivariate_gaussian_pdf(x, mu, sigma, epsilon=1e-6):
    D = len(mu)
    x = np.asarray(x)
    
    # 添加小常数防止奇异矩阵
    sigma_reg = sigma + epsilon * np.eye(D)
    
    # 计算行列式和逆（使用伪逆更安全）
    det = np.linalg.det(sigma_reg)
    inv_sigma = np.linalg.pinv(sigma_reg)
    
    # 对数空间计算更稳定
    log_norm_const = -0.5 * (D * np.log(2*np.pi) + np.log(det))
    x_mu = x - mu
    log_exponent = -0.5 * x_mu.T @ inv_sigma @ x_mu
    
    return np.exp(log_norm_const + log_exponent)

3. 计算区域概率：从理论到实践

计算单个点的概率密度相对简单，但实际应用中，我们更常需要计算随机变量落在某个区域的概率。对于多元高斯分布，这需要计算多重积分：

$$ P(a \leq x \leq b) = \int_{a_1}^{b_1} \int_{a_2}^{b_2} f_{\mu,\Sigma}(x) dx_2 dx_1 $$

解析解通常难以求得，因此我们转向数值方法。以下是几种实用方法：

3.1 蒙特卡洛积分法

蒙特卡洛方法通过随机采样来估计积分值。基本步骤是：

1. 在感兴趣区域内生成大量随机点
2. 计算这些点处的概率密度
3. 取平均值并乘以区域体积

def monte_carlo_integration(mu, sigma, lower_bounds, upper_bounds, n_samples=100000):
    D = len(mu)
    
    # 1. 在矩形区域内生成均匀随机样本
    samples = np.random.uniform(low=lower_bounds, high=upper_bounds, 
                               size=(n_samples, D))
    
    # 2. 计算每个样本的概率密度
    pdf_values = np.array([multivariate_gaussian_pdf(x, mu, sigma) 
                          for x in samples])
    
    # 3. 计算区域体积
    volume = np.prod(upper_bounds - lower_bounds)
    
    # 4. 估计积分值
    integral_estimate = volume * np.mean(pdf_values)
    
    return integral_estimate

让我们计算身高[70,80]英寸，体重[65,75]公斤范围内的概率：

lower = np.array([70, 65])  # 身高下限，体重下限
upper = np.array([80, 75])  # 身高上限，体重上限

prob = monte_carlo_integration(mu, sigma, lower, upper)
print(f"Probability using Monte Carlo: {prob:.4f}")

3.2 网格近似法

对于低维问题（如我们的二维案例），网格法可能更高效：

def grid_integration(mu, sigma, lower_bounds, upper_bounds, step=0.1):
    # 创建网格
    x1 = np.arange(lower_bounds[0], upper_bounds[0], step)
    x2 = np.arange(lower_bounds[1], upper_bounds[1], step)
    xx1, xx2 = np.meshgrid(x1, x2)
    
    # 计算每个网格点的概率密度
    positions = np.vstack([xx1.ravel(), xx2.ravel()]).T
    pdf_values = np.array([multivariate_gaussian_pdf(x, mu, sigma) 
                          for x in positions])
    
    # 计算积分近似值
    integral = np.sum(pdf_values) * (step ** len(mu))
    
    return integral

比较两种方法的结果：

print(f"Probability using grid method: {grid_integration(mu, sigma, lower, upper):.4f}")

3.3 使用SciPy的数值积分

对于更精确的计算，可以使用SciPy的数值积分功能：

from scipy.integrate import nquad

def integrand(x1, x2):
    x = np.array([x1, x2])
    return multivariate_gaussian_pdf(x, mu, sigma)

# 定义积分界限
bounds = [(70, 80), (65, 75)]

# 计算数值积分
result, error = nquad(integrand, bounds)
print(f"Scipy numerical integration: {result:.4f} (error: {error:.2e})")

4. 实际应用中的注意事项与优化

在实际项目中应用这些技术时，有几个关键点需要注意：

4.1 协方差矩阵的正定性

协方差矩阵必须是正定的。在实践中，我们可能会遇到：

• 样本量不足导致的秩缺失
• 数值误差造成的微小负特征值

解决方法包括：

def make_positive_definite(sigma, epsilon=1e-6):
    """
    通过添加小常数到对角线确保矩阵正定
    """
    n = sigma.shape[0]
    return sigma + epsilon * np.eye(n)

4.2 高维诅咒

随着维度增加，概率质量越来越集中在边缘区域，导致数值积分变得困难。这时可以考虑：

• 重要性采样
• 马尔可夫链蒙特卡洛(MCMC)方法
• 变分近似

4.3 对数概率计算

对于非常小的概率值，直接计算可能导致数值下溢。改用对数空间更稳定：

def log_multivariate_gaussian_pdf(x, mu, sigma):
    D = len(mu)
    x = np.asarray(x)
    x_mu = x - mu
    inv_sigma = np.linalg.inv(sigma)
    
    log_det = np.log(np.linalg.det(sigma))
    log_norm = -0.5 * (D * np.log(2*np.pi) + log_det)
    log_exp = -0.5 * x_mu.T @ inv_sigma @ x_mu
    
    return log_norm + log_exp

4.4 批量计算优化

当需要计算大量点的概率密度时，使用向量化操作可以大幅提升性能：

def batch_multivariate_gaussian_pdf(X, mu, sigma):
    """
    X: (N,D)矩阵，N个D维样本
    mu: (D,)均值向量
    sigma: (D,D)协方差矩阵
    """
    D = mu.shape[0]
    X = np.asarray(X)
    inv_sigma = np.linalg.inv(sigma)
    log_det = np.log(np.linalg.det(sigma))
    
    X_mu = X - mu[None,:]
    exponent = -0.5 * np.sum(X_mu @ inv_sigma * X_mu, axis=1)
    log_norm = -0.5 * (D * np.log(2*np.pi) + log_det)
    
    return np.exp(log_norm + exponent)

5. 完整案例：身高体重概率计算

现在，让我们将所有内容整合到一个完整的案例中。假设我们有以下任务：

给定身高体重服从N(μ,Σ)分布，其中： μ = [75.0, 71.3] Σ = [[874, 327], [327, 929]] 计算身高在[70,80]英寸且体重在[65,75]公斤范围内的概率

5.1 准备工作

首先，我们实现一个综合解决方案：

import numpy as np
from scipy.stats import multivariate_normal
from scipy.integrate import nquad

class MultivariateGaussianCalculator:
    def __init__(self, mu, sigma):
        self.mu = np.asarray(mu)
        self.sigma = np.asarray(sigma)
        self.D = len(mu)
        self.validate_parameters()
        
    def validate_parameters(self):
        assert self.sigma.shape == (self.D, self.D), "Sigma must be DxD matrix"
        assert np.allclose(self.sigma, self.sigma.T), "Sigma must be symmetric"
        
        # 确保正定
        try:
            np.linalg.cholesky(self.sigma)
        except np.linalg.LinAlgError:
            self.sigma = self._regularize_sigma()
    
    def _regularize_sigma(self, epsilon=1e-6):
        return self.sigma + epsilon * np.eye(self.D)
    
    def pdf(self, x):
        """计算单点概率密度"""
        x = np.asarray(x)
        inv_sigma = np.linalg.inv(self.sigma)
        x_mu = x - self.mu
        exponent = -0.5 * x_mu.T @ inv_sigma @ x_mu
        norm = 1.0 / ((2 * np.pi) ** (self.D/2) * np.sqrt(np.linalg.det(self.sigma)))
        return norm * np.exp(exponent)
    
    def region_probability(self, lower_bounds, upper_bounds, method='scipy', **kwargs):
        """
        计算区域概率
        method: 'scipy', 'monte_carlo', 或 'grid'
        """
        methods = {
            'scipy': self._scipy_integrate,
            'monte_carlo': self._monte_carlo_integrate,
            'grid': self._grid_integrate
        }
        return methods[method](lower_bounds, upper_bounds, **kwargs)
    
    def _scipy_integrate(self, lower_bounds, upper_bounds):
        def integrand(*args):
            x = np.array(args)
            return self.pdf(x)
        
        bounds = list(zip(lower_bounds, upper_bounds))
        result, _ = nquad(integrand, bounds)
        return result
    
    def _monte_carlo_integrate(self, lower_bounds, upper_bounds, n_samples=100000):
        samples = np.random.uniform(low=lower_bounds, high=upper_bounds,
                                  size=(n_samples, self.D))
        pdf_values = np.apply_along_axis(self.pdf, 1, samples)
        volume = np.prod(np.array(upper_bounds) - np.array(lower_bounds))
        return volume * np.mean(pdf_values)
    
    def _grid_integrate(self, lower_bounds, upper_bounds, step=0.1):
        axes = [np.arange(l, u, step) for l, u in zip(lower_bounds, upper_bounds)]
        grid = np.meshgrid(*axes)
        positions = np.vstack([g.ravel() for g in grid]).T
        pdf_values = np.apply_along_axis(self.pdf, 1, positions)
        return np.sum(pdf_values) * (step ** self.D)

5.2 计算结果比较

现在，我们可以用不同方法计算并比较结果：

# 初始化计算器
mu = [75.0, 71.3]
sigma = [[874, 327], [327, 929]]
calculator = MultivariateGaussianCalculator(mu, sigma)

# 定义积分区域
lower = [70, 65]
upper = [80, 75]

# 使用不同方法计算
methods = ['scipy', 'monte_carlo', 'grid']
results = {}

for method in methods:
    prob = calculator.region_probability(lower, upper, method=method)
    results[method] = prob
    print(f"{method:12s}: {prob:.6f}")

# 可视化比较
plt.bar(results.keys(), results.values())
plt.ylabel('Probability')
plt.title('Comparison of Integration Methods')
for method, prob in results.items():
    plt.text(method, prob, f"{prob:.4f}", ha='center', va='bottom')
plt.show()

5.3 结果分析与选择

三种方法得到的结果应该非常接近，但各有优缺点：

方法	优点	缺点	适用场景
SciPy积分	最精确，误差可控	计算量大，高维时不可行	低维精确计算
蒙特卡洛	实现简单，维度扩展性好	结果有随机性，需要大量样本	中高维近似计算
网格法	确定性结果，实现直观	维度灾难，高维时计算量爆炸	超低维(2-3维)快速计算

在实际项目中，选择哪种方法取决于具体需求。对于我们的身高体重案例（2维），三种方法都适用，SciPy积分可能是最佳选择。