高斯分布概率密度的二重积分：从查表到真正算出置信椭圆的完整过程
我第一次认真面对“高斯分布概率密度的二重积分”这个题目，是在给一批二维定位数据画置信椭圆的时候。当时查了不少资料，一维正态分布有现成的经验法则，均值加减三个标准差对应99.73%的概率，这个谁都会说。可一旦到了二维，问题就变成了“椭圆内部的概率到底是多少”，市面上很多文章直接给一张表格了事，从来没说这张表是怎么来的，更没人告诉你当椭圆倾斜、变量相关、根本不是一个规整圆的时候，这张表还怎么用。后来被现实反复教育过才明白，要回答这些问题，至少得把二维高斯密度的二重积分从头到尾想清楚。
这篇文章写的是我在这件事上的完整记录：先从一维到二维的推广讲起，点出二重积分到底卡在哪里；然后把最实用的计算路径拆开，包括极坐标变换、标准化、马氏距离、卡方分布这条捷径；最后给出能直接跑的数值算法和代码，再把高维、相关变量的坑一并讲清楚。适合做数据分析、信号处理、定位导航、传感融合或者统计建模的朋友参考，尤其是那些画过置信椭圆、出过误差分布图、却还没搞明白“椭圆里那一坨概率到底怎么算”的人。
1. 从一维查表到二维查不了表：二重积分到底卡在哪里
1.1 一维的时候你只需要一个数，二维的时候你需要一整张曲面
一维高斯分布的概率密度函数大家都熟：
f(x) = (1 / (σ√(2π))) · exp(−(x−μ)² / (2σ²))
求某个区间 [a, b] 内的概率，就是对这个密度函数做一维积分：
P(a ≤ X ≤ b) = ∫[a,b] f(x) dx
这个积分没有初等原函数，但数值上已经解决得透透的。学过统计的都知道，标准化之后查标准正态分布表，或者直接调用统计软件里的正态累积分布函数（CDF），一秒钟出结果。查表查了这么多年，从没觉得有什么不对劲——直到变量从一个变成两个。
二维高斯分布的概率密度不再是曲线下面的面积，而是曲面下面的“体积”。联合密度函数写成：
f(x, y) = 1 / (2π |Σ|^{1/2}) · exp(−½ (z−μ)ᵀ Σ⁻¹ (z−μ))
其中 z = (x, y)ᵀ，Σ 是二维协方差矩阵，|Σ| 是行列式。要求“数据落入某个区域 D 的概率”，就是算：
P((X, Y) ∈ D) = ∬_D f(x, y) dx dy
麻烦就麻烦在这里：一维的 CDF 有一个成熟的函数可以调用，但二维的“联合 CDF”没有一个简单的表达式可以解析写出，尤其当关注的是一个椭圆或者任意形状区域时，几乎所有情况下都得老老实实做二重积分。
我遇到过很多同事，印象里还停留在一维那套“均值±3σ”的思维。到了二维，上来就问：“两个维度各自在 ±3σ 范围内，概率是不是也约等于99.73%？”答案是：只有两个维度完全独立、并且你恰好关心的是一个与坐标轴对齐的矩形区域时，才能用一维概率相乘；一旦关心的是椭圆区域，或者变量之间有相关性，那个“69.7%×3σ”的说法就根本不成立了。下面会仔细算给你看。
1.2 工程里最常见的二维积分问题：置信椭圆与圆概率误差
二维高斯积分在工程里出现频率最高的地方，就是置信椭圆，以及定位领域里经常提到的圆概率误差（CEP，Circular Error Probable）。
先说置信椭圆。任何一组服从多元正态分布的二维数据，都可以利用协方差矩阵的特征值和特征向量画一个椭圆。椭圆的中心放在均值点，长轴方向由最大特征值对应的特征向量决定，长轴半径取特征值的平方根乘以某个缩放系数。这个缩放系数就是你想要的置信水平。比如你要一个“95%置信椭圆”，就需要找到一个缩放系数 t，使得：
P(马氏距离 ≤ t) = 0.95
这里的概率计算本质上就是高斯密度在一个椭圆区域上的二重积分。椭圆的形状由协方差矩阵决定，缩放系数是积分区域的上限，积分密度是联合高斯密度。两者一组合，面积覆盖和概率覆盖就对应上了。
再说圆概率误差。在定位系统评估中，经常要回答“误差有一半概率落在半径多大的圆里”这个问题。假设定位误差在两个正交方向上独立同分布，都是零均值、标准差 σ 的正态分布，那么误差向量 (x, y) 的联合密度就是两个一维密度的乘积。要求“落进半径为 r 的圆的概率达到50%”，就是求：
∬_{x²+y² ≤ r²} (1/(2πσ²)) exp(−(x²+y²)/(2σ²)) dx dy = 0.5
解出这个 r，然后把 r 除以 σ 得到比例系数，就是定位领域里常见的CEP系数。这个系数在不同导航、雷达评估里都有固定值，但值的来历就是上面这个二重积分。我最早见到这个题目的场景，恰好就是在给一组测量设备做误差半径评估，那时候还不理解为什么概率50%对应的圆半径约等于1.177倍的σ。后来自己动手把这个积分拆开算了一遍，才彻底明白那个1.177是怎么冒出来的。
2. 二重积分的核心直觉：极坐标变换把“二维查表”变成“一维查表”
2.1 为什么高斯密度在极坐标下会突然变“好算”
一个让很多人感到神奇的细节：一维高斯密度积分没有初等表达式，但二维高斯密度在全平面上的积分——也就是 (1/(2π)) ∫∫ exp(−(x²+y²)/2) dx dy ——却可以精确等于1，而且证明过程只需要几步极坐标变换。
具体做法是这样的。把直角坐标换成极坐标：
x = r cosθ, y = r sinθ
面积元 dx dy 变成 r dr dθ。于是：
∫∫ exp(−(x²+y²)/2) dx dy = ∫₀^{2π} ∫₀^{∞} exp(−r²/2) · r dr dθ
内层对 r 的积分，正好是 ∫ exp(−r²/2) r dr = 1，外层 2π 一乘，再除以归一化常数 2π，结果就是1。
这个“r 恰好出现”不是巧合。核心原因是二维高斯密度的指数部分只依赖半径 r，也就是距离原点的长度，不依赖方向；而极坐标的面积元恰好带一个 r 的权重因子，它和指数函数的导数形状配合得天衣无缝。θ 方向是均匀的，可以直接积掉；只剩 r 方向一个一维积分。这个性质几乎贯穿了所有二维高斯积分的快速算法。
这个结论在直观上也很重要：标准二维高斯分布（零均值、单位方差、相互独立）在任何以原点为圆心的圆形区域内，落进去的概率只和半径有关。设半径为 R，概率就是：
P(r ≤ R) = 1 − exp(−R²/2)
这个公式简直可以用来当面试题。R=1 时概率约0.3935，R=√2≈1.414 时概率约0.6321，R=2 时约0.8647，R=3 时约0.9889。注意，二维跑到R=3（即两个方向各偏离3σ）时，概率是98.89%，不是一维的99.73%。这个差值在工程上是实打实的，画3σ圆却想标99.73%置信度，是小细节，但错了就是错了。
2.2 从圆域到椭圆域：“缩放”和“旋转”不改变问题的本质
如果有协方差矩阵，问题就从“圆形区域内的积分”变成“椭圆区域内的积分”了。第一反应是头大，但线性代数帮了大忙。
高斯密度的指数部分是马氏距离的平方：
D² = (z−μ)ᵀ Σ⁻¹ (z−μ)
协方差矩阵 Σ 是实对称正定矩阵，可以做特征值分解 Σ = UΛUᵀ，U 是正交矩阵（旋转），Λ 是对角矩阵（缩放）。做一个坐标变换：
w = Λ^(−1/2) Uᵀ (z−μ)
这个变换把原坐标系下的椭圆变成了一个单位圆。因为 U 对应旋转，Λ^(−1/2) 对应把各轴按特征值的倒数缩放到单位尺度，所以：
- 原空间里的置信椭圆 z 满足 D² ≤ t²；
- 新空间里的点 w 满足 |w| ≤ t，也就是一个半径为 t 的圆。
而且雅可比行列式算下来，恰好把密度函数里的 1/(2π |Σ|^{1/2}) 抵消掉，新空间里就是一个标准二维高斯分布。换句话说：任何二维高斯分布的二重积分，都可以通过“旋转+缩放”变成标准二维高斯在圆形区域上的积分。
这就是我文章标题里那个“二重积分”的真正解法思路。它没有直接用数值方法硬刚椭圆区域，而是先在纸上把问题化简成标准形状，再用上一步的极坐标结果。
实际操作里，我从不建议直接把椭圆参数填进积分器里硬算。正确姿势是：先标准化，再积分，最后用卡方分布出概率。下面第三节详细拆这一步。
3. 先标准化，再积分：把一般二重积分变成马氏距离的查表问题
3.1 协方差矩阵不是零，你的“积分区域”就必须重新画
很多从一维思维转到二维的人，最喜欢直接把积分区域写成“x 从 μ±aσ_x、y 从 μ±bσ_y”的矩形。这个做法在变量独立时是对的，一旦变量之间有相关性，比如协方差矩阵的对角线之外有非零项，问题就来了：椭圆是倾斜的，但矩形和坐标轴对齐，矩形边界既不是椭圆切线，又不能反映相关方向的信息。
举个例子。假设 Σ = [[1, 0.7], [0.7, 1]]，强正相关，均值在原点。你画一个矩形 [-1, 1] × [-1, 1]，并在矩形里积分，得到的概率大约是0.2864。可如果你真正想要的是“1个标准差级别的椭圆”概率，真实值不是这个，而是约0.3935（因为标准圆半径1时概率为 1 − exp(−0.5)）。这两个数差得远，原因就在于矩形把椭圆尖角区域漏掉、又把椭圆外沿区域包进来了，一加一减不是互相抵消，通常偏小。
所以做二维高斯积分，第一步永远是：把积分区域从“直角坐标矩形”翻译成“马氏距离约束”。凡是能用马氏距离 D² ≤ c 描述的椭圆区域，就被这个约束唯一决定了。做数值积分之前先把 c 算对，比考虑用什么积分方法都重要。
3.2 椭圆区域内概率：用卡方分布一招搞定
一旦确定了马氏距离阈值 c，概率可以直接用自由度为2的卡方分布计算。因为新空间里 w₁、w₂ 是两个独立标准正态变量，w₁² + w₂² 恰好服从自由度2的卡方分布。然后：
P(D² ≤ c) = P(χ²₂ ≤ c) = 1 − exp(−c/2)
以前面那个“置信椭圆”为例：要画一个覆盖概率 p 的椭圆，缩放系数 t 就是：
t = sqrt(χ²_{2,p})
95% 对应 χ²₂,0.95 ≈ 5.9915，所以 t ≈ 2.4477；99% 对应 χ²₂,0.99 ≈ 9.2103，t ≈ 3.0349。这个 t 就是标准单位圆半径，再还原到原坐标系就要乘上对应的特征值平方根。
第一次算的时候我有点怀疑：极坐标下那个 R 和卡方分布的 d 是不是一回事？后来想明白了，极坐标里 r² = x²+y²，对标准二维高斯就是卡方分布的自由度2随机变量。所以极坐标公式 1−exp(−R²/2) 只是卡方分布CDF在具体R值处的写法，两者完全等价。查表时代用卡方分布表，手算时代用指数公式，本质是同一个数学事实。
3.3 矩形区域为什么没有捷径：只能回到数值积分
如果积分区域是任意矩形，那就没有卡方分布这种捷径了。矩形边界把 θ 方向上的均匀性破坏了，极坐标变换后 θ 的积分范围不再是完整一周，r 的上限也随 θ 变化，成了一个分段函数。这种带边界的积分没有统一的封闭解，只能上数值方法。
最常见的场景是“双变量CDF”。比如要算 P(X ≤ a, Y ≤ b)，积分区域是 (-∞, a] × (-∞, b]，这就是一个右上角半无限矩形。它的数值积分有现成算法，大多数统计工具里都有“二元正态累计概率”这种函数。但很多需求并不是标准矩形，而是“x 在 [−1, 1]、y 在 [0.5, 2]”这种有限矩形。这类问题没有现成的一行函数，老老实实做二重数值积分是最稳的方案。
我自己的经验是：先判断积分区域形状。如果是椭圆/球形区域，一律先标准化，再用卡方分布算，又快又准；如果是矩形或不规则区域，再考虑数值积分。这个判断顺序能帮你省下大量时间，也避免把椭圆积分硬塞给数值积分器然后发现精度不够的尴尬。
4. 数值求解实操：用Python把二重积分真正跑起来
4.1 三个可选方法的威力边界：网格法、蒙特卡洛法、自适应积分法
当问题确实需要数值二重积分时，我常在这三类方法里选：
第一类是等距网格矩形法。把积分区域划分成小格，每个小格取中心点，用密度值乘以小格面积求和。实现简单，能出“大概”的结果，但二维网格的主要问题是精度和计算量的矛盾。要把误差压到千分之一量级，网格密度得非常高，在二维还能忍受，三维以上就崩盘了。
第二类是蒙特卡洛法。在积分区域内随机撒 N 个点，统计落在区域内的比例乘以区域面积，或者直接对密度函数在采样点上取均值。它的优势是维数高了照样能跑，误差随 N^(−1/2) 下降，从千分之一精度到万分之五精度大概需要4倍点数。劣势是结果带随机性，做验收时不稳定，同一段代码跑两次可能差出小数点后第三位。
第三类是自适应积分法，也是我遇到困难时的首选。像数值积分工具里的 dblquad 这类函数，底层是自适应高斯-勒让德或者类似算法，会按积分的局部行为自动加密采样点，精度高且稳定。代价是代码稍复杂、当被积函数非常尖或者积分区域边界复杂时可能会慢，但二维场景完全可接受。
对二维高斯密度这种光滑、快速衰减的函数，我给的建议很明确：能用解析法/卡方法就先用；需要数值时优先自适应积分；蒙特卡洛作为验证手段或高维兜底。
4.2 代码示例：标准二维高斯在矩形区域的积分
下面用一个实际例子跑一遍。问题：计算标准二维高斯分布（零均值、单位方差、独立）在单位矩形 [0,1]×[0,1] 内的概率。
import numpy as np from scipy.integrate import dblquad from scipy.stats import norm def std_gaussian_2d(y, x): return (1.0 / (2.0 * np.pi)) * np.exp(-0.5 * (x*x + y*y)) # 区域：x in [0,1], y in [0,1] prob, err = dblquad( std_gaussian_2d, 0.0, 1.0, # x 下限、上限 lambda x: 0.0, lambda x: 1.0 # y 下限、上限 ) # 解析对照：两个独立一维正态分别在 [0,1] 内的概率相乘 analytic = (norm.cdf(1.0) - norm.cdf(0.0))**2 print(f"数值积分: {prob:.10f}, 误差估计: {err:.2e}") print(f"一维相乘: {analytic:.10f}")因为两个维度独立，这里可以用一维CDF乘积做对照。运行结果数值积分和一维相乘应当一致到小数点后八九位。如果你的维度之间不独立，就不能这么对照了，那就需要下面第4.3节的标准化处理。
4.3 代码示例：有相关性的高斯分布怎么算任意区域积分
有相关性的情况，核心是先对协方差矩阵做Cholesky分解或者特征分解，把变量变成“标准高斯组合”，再去数值积分。
import numpy as np from scipy.integrate import dblquad from scipy.stats import chi2 mu = np.array([0.0, 0.0]) Sigma = np.array([[1.0, 0.7], [0.7, 1.0]]) L = np.linalg.cholesky(Sigma) def standard_gaussian_pdf(w1, w2): return (1.0 / (2.0 * np.pi)) * np.exp(-0.5 * (w1*w1 + w2*w2)) # 在 w 空间（标准高斯空间）定义原矩形 # 原矩形假设: x in [-1,1], y in [-1,1] # Jacobian: dx dy = |det L| dw1 dw2 = sqrt(det Sigma) det_jac = np.linalg.det(L) def rect01(w1, w2, a=1.0): # 原坐标系下的点 z = mu + L @ np.array([w1, w2]) x, y = z[0], z[1] # 判断是否在矩形内 return 1.0 if (-a <= x <= a and -a <= y <= a) else 0.0 # 直接在标准高斯空间积分 prob, err = dblquad( lambda w2, w1: standard_gaussian_pdf(w1, w2) * rect01(w1, w2), -8.0, 8.0, lambda w1: -8.0, lambda w1: 8.0 ) prob = prob * det_jac print(f"矩形[-1,1]x[-1,1]内的概率: {prob:.6f}") # 对照：马氏距离阈值 c=1 的椭圆概率 c = 1.0 prob_ellipse = chi2.cdf(c, df=2) print(f"马氏距离<=1（即单位椭圆）: {prob_ellipse:.6f}")代码里有个细节值得解释：我在标准高斯空间里用“指示函数”判断每个点是否落在原矩形内，再乘上雅可比行列式。这样就不需要特意推导被积函数在直角坐标下的表达式，直接把区域边界翻译成判断条件就行。积分上下限取正负8，是因为标准高斯在±8之外的密度已经完全可忽略，不会对积分结果产生可观察影响；如果你对精度抠得比较细，可以把边界再拉宽到±10，代价就是多算几轮，影响不大。
跑完之后你会发现：矩形 [-1,1]×[-1,1] 里的概率大约是0.2864，而马氏距离≤1的椭圆概率是0.3935。这两个值差了不少，正好回到第3.1节那个提醒——区域形状不一样，概率完全不同，别混为一谈。
5. 高维推广与几个我替你先踩过的坑
5.1 维度升高之后：网格法先完蛋，蒙特卡洛续命
高斯密度积分从二维往高维推广，思路是相同的：标准化之后，任何多元高斯分布在马氏距离球形区域内的概率都可以用自由度为 d 的卡方分布计算，P(D² ≤ c) = χ²_d(c)。这是高维场景下最快的路。
但如果你非要做数值积分，维数升高意味着网格法马上失灵。三维网格每个维度分100份，就是100万个点，虽然也能跑，但到四维五维就指数爆炸了。这时候蒙特卡洛的相对优势显现出来：误差和维数关系不大，只和采样量有关。五维十维的场景，我一般先用卡方分布算球形概率，形状不规则的区域再用重要性采样或者直接蒙特卡洛硬撒点。撒点数量到十万一百万量级，统计上基本够用。
5.2 坑一：以为“两个维度独立”就可以直接相乘，结果算漏了相关性
这是我在实际项目里看到最多、也最容易犯的错。一组数据做二维高斯拟合，协方差矩阵的非对角项明明不是零，某同学把二维CDF简写成两个一维CDF的乘积。相关性为0.7时，这样算出来的概率误差非常明显，区间概率算小或者算大，取决于区域位置和符号。
为什么不能简单地相乘？因为两个相关正态变量的联合密度不是一个乘积形函数。它必须包含协方差项，表达式里有一个 exp(−(x²−2ρxy+y²)/(2(1−ρ²))) 形状的因子。相关性越强，联合密度在一条对角线附近越集中，矩形区域里的概率和一维乘积的差异越大。要安全地处理，要么用之前说的标准化+指示函数方法，要么用协方差矩阵直接构造联合密度再数值积分，反正不能想当然地拆成两个独立变量。
5.3 坑二：蒙特卡洛结果的抖动被当成算法不稳定
蒙特卡洛跑完一次出结果，再跑一次差零点零零几，有同事以为代码写错了。这不是错误，是蒙特卡洛的正常方差。验收时如果想看“确定性结果”，应该用自适应积分法再算一次，两个结果在误差范围内对齐就够了。蒙特卡洛的价值在于快速试探、高维兜底，不适合当作唯一答案输出。尤其在做置信椭圆这种需要精确边界的任务里，我建议以自适应积分结果为准，蒙特卡洛只做交叉验证。
5.4 坑三：积分区域边界设得太窄，结果悄悄偏小
有次帮人排查一个二维高斯积分，对方把积分上下限设在正负3，以为3σ之外概率可以忽略。在“全平面都能积”的一维思维里，正负3确实覆盖了99.73%——但注意，那是单变量从负无穷积到正无穷、然后取区间内部的结果。你要是真的把积分边界截断在±3再去算一定区域的概率，截断误差是累计在结果里的。对矩形区域还好，正负3的边界截断对标准高斯带来的误差约在10的负几次方量级；但对高维或者偏态协方差矩阵，尾部的重要性会放大，最好把积分的“外部背景区域”放宽到±6甚至更高，保证数值积分器在真正“算完”整个有效质量时不会漏掉尾部。
6. 工程落地时要带走的几个结论
二重积分本身是数学工具，但在工程里它永远连着具体参数。数据可视化时要画置信椭圆，你得先定概率水平，再反推缩放系数；定位评估里要报CEP，你得先做二重积分解出半径系数；滤波算法里算更新门限，本质也是一个高斯分布球域积分问题。这些场景里，我已经习惯了一套固定的处理顺序：先判断区域形状，能标准化就标准化，能用卡方分布就用卡方分布，实在不规则才上数值积分。
有些细节值得反复说。第一，标准化变换矩阵要保留下来，不要算完概率就丢掉。因为很多时候你不光要概率，还要把标准空间的圆还原成原空间的椭圆去画图。有了协方差矩阵的特征值分解，旋转和缩放都随时能恢复，省得二次开发时重新推一遍。第二，做验收前先跑一个简单解析对照。比如标准二维高斯在单位圆内的概率是 1−exp(−0.5)，在 [-1,1]×[-1,1] 矩形内是 (erf(1/√2))²，这些现成结果能快速验证你的代码有没有基本错误。第三，一定要把“矩形区域”和“椭圆区域”分开回答。很多需求文档里说“给我二维高斯95%置信区间”，实际上他们想要的是置信椭圆，这时候你用矩形积分算出来的数字填进去，可能在可视化时看起来是对的，但概率含义完全不对。
我个人在实际项目里的体会是：这个题目最容易卡人的不是积分本身，而是没有先意识到二维和一维的差异。一维的“三西格玛法则”太顺口了，顺口到让人下意识在二维里也按同样的值报概率。真正动手做一遍极坐标变换、标准化、卡方分布查表或者数值积分之后，你不仅会算置信椭圆了，更会对“相关”“协方差”“区域形状”这些概念产生一种新的直觉——知道哪里该取捷径，哪里必须老老实实上数值工具。这也是我写这篇文章最主要的动机：把这条从公式到代码、从标准圆到椭圆再到任意区域的路完整走一遍，让后来者少花几个晚上在摸索上。
