高斯积分这个名字很多人是在大学物理课上第一次听到的但真正被它咬一口往往是在写代码或者对公式的时候。前阵子一个做分子运动模拟的同学拿着速率分布曲线来找我说怎么画都归一化不到 1公式抄了三遍还是不对最后发现是 $\int_0^\infty x^2 e^{-ax^2},dx$ 那个系数里的 $a^{3/2}$ 被写成了 $a^{1/2}$。这种错误很典型不是不会是没真正把高斯积分这一族式子当成一个整体去理解只记住了最常用的那一个 $\sqrt{\pi}$。高斯积分本质上是一类以 $e^{-x^2}$ 为核心的反常积分从形式上看它平平无奇但它是高等数学里反常积分技巧的集大成者是大学物理里麦克斯韦速率分布、量子谐振子波函数、热传导核、配分函数的标准工具也是概率统计里正态分布归一化常数、特征函数、多维高斯密度、高斯过程的底层来源。换句话说凡是出现钟形曲线的地方背后都站着一个已经算好的高斯积分。这篇文章面向三类人正在学高数但被反常积分绕晕的低年级学生、做物理或统计建模需要反复调用这些公式的研究者、以及写代码时需要做数值归一化和概率计算但又不想每次都从头推一遍的工程实现者。我不打算把它写成公式大全那样网上到处都是。我更想讲清楚的是三件事这个积分为什么非得绕个弯才能算出来它的参数化家族为什么才是真正能用的形态以及在代码里实现它的时候那些教材不会提的坑。1. 先搞清楚它到底难在哪为什么这个积分不能用牛顿-莱布尼茨公式1.1 找不到初等原函数不代表积不出来初学者最容易产生的困惑是既然 $\int_0^1 x^2,dx$ 能算为什么 $\int_0^\infty e^{-x^2},dx$ 就算不出来答案藏在原函数的存在性里。$e^{-x^2}$ 在任意有限区间上当然连续、当然可积但它的原函数没法用有限个初等函数多项式、指数、对数、三角、反三角及其复合表示出来。这一点是刘维尔在十九世纪就证明过的结论不是没人聪明是数学上不可能。所以在牛顿-莱布尼茨公式这条路上你会被卡死在第一步找不到 $F(x)$ 就代不了上下限。工程上常用的做法是定义一个新函数来命名这个原函数也就是误差函数 $\operatorname{erf}(x) \frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2},dt$然后把它当作查表对象。但注意这只是把问题改了名字并没有真正解决。真正的突破点是虽然原函数写不出来可它在 $(-\infty, \infty)$ 上的定积分却有一个极其漂亮的闭式结果也就是 $\sqrt{\pi}$。这里有个直觉值得建立起来。$e^{-x^2}$ 衰减得比任何 $e^{-|x|}$ 都快但又比 $e^{-x^4}$ 慢这个中间速度恰好和圆周率挂上了钩。为什么会冒出 $\pi$一个本来只跟一维直线有关的积分结果里却有圆的性质这提示我们解法一定和一个二维的、带旋转对称性的构造有关。提示不要在一个方向里死磕这类积分的解法都建立在升维上一维走不通就跳到二维。1.2 平方一下一维难题降维成二维送分题这是整个推导里最漂亮的一步。设 $I \int_{-\infty}^{\infty} e^{-x^2},dx$。因为被积函数恒正所以 $I 0$且积分收敛因此可以平方$$I^2 \left(\int_{-\infty}^{\infty} e^{-x^2},dx\right)\left(\int_{-\infty}^{\infty} e^{-y^2},dy\right) \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} e^{-(x^2y^2)},dx,dy$$这一步的合法性来自富比尼定理两个积分都绝对收敛所以可以写成二重积分并且两个一维积分的结果可以相乘。这一步很多人会跳过觉得显然但严格性恰恰在这里。接下来换成极坐标$x r\cos\theta$$y r\sin\theta$$x^2 y^2 r^2$而面积元 $dx,dy$ 要换成 $r,dr,d\theta$其中那个多出来的 $r$ 是雅可比行列式的值。这是极坐标换元的固定规则很多人抄公式时会忘记 $r$然后就得出一个错误的结果。换元后$$I^2 \int_0^{2\pi}\int_0^\infty e^{-r^2} r,dr,d\theta 2\pi \cdot \left[-\frac{1}{2}e^{-r^2}\right]_0^\infty 2\pi \cdot \frac{1}{2} \pi$$于是 $I \sqrt{\pi}$。整个过程的关键在于$e^{-(x^2y^2)}$ 这个函数具有旋转对称性在直角坐标下两个变量纠缠不清在极坐标下直接拆成了两个独立的部分内层关于 $r$ 的积分因为多了那个 $r$ 因子正好可以用换元 $u r^2$ 一步做出来。如果少了极坐标里那个 $r$看起来形式一样但内层积分就变成了 $\int e^{-r^2}dr$又回到原点。注意这个技巧只对整个实数轴成立。如果换成 $\int_0^\infty$用对称性可以直接得到半个结果 $\sqrt{\pi}/2$但如果上限是有限数比如 $a$那就没有闭式解了只能用误差函数表示。我个人的体会是带学生推这一步时最值钱的地方不是结果而是让他们自己意识到什么时候该考虑极坐标。判据很简单被积函数里出现了 $x^2y^2$ 这种形式的整体或者积分区域是圆盘、全平面那就该往极坐标想。1.3 记住一个结果没用记住一族结果才够用只记住 $\int_{-\infty}^{\infty}e^{-x^2}dx \sqrt{\pi}$在真正做题时是远远不够的。实际遇到的高斯积分几乎从来不是标准形式一定带着各种参数有 $\sigma$、有 $a$、有 $m/(2kT)$、有 $1/(4Dt)$。所以真正需要内化的是一个参数化家族核心是下面这几个积分形式结果常见用途$\int_{-\infty}^{\infty} e^{-ax^2}dx$$\sqrt{\pi/a}$归一化常数的基础形$\int_0^{\infty} e^{-ax^2}dx$$\frac{1}{2}\sqrt{\pi/a}$半轴上的分布$\int_{-\infty}^{\infty} e^{-ax^2bx}dx$$\sqrt{\pi/a},e^{b^2/(4a)}$特征函数、配方法$\int_0^{\infty} x,e^{-ax^2}dx$$\frac{1}{2a}$极坐标径向积分、光束功率$\int_0^{\infty} x^2e^{-ax^2}dx$$\frac{\sqrt{\pi}}{4a^{3/2}}$速率分布、方差$\int_0^{\infty} x^3e^{-ax^2}dx$$\frac{1}{2a^2}$平均速率$\int_0^{\infty} x^4e^{-ax^2}dx$$\frac{3\sqrt{\pi}}{8a^{5/2}}$方均根、四阶矩这张表里的所有 $a$ 都要求 $a 0$否则积分发散这一点绝对不能漏。检查方法很直接如果 $a \le 0$$e^{-ax^2}$ 在无穷远处不趋于零甚至增大积分必然发散。统一记法是把它们全部写成 Gamma 函数的形式$$\int_0^\infty x^{n}e^{-ax^2},dx \frac{1}{2},\frac{\Gamma!\left(\frac{n1}{2}\right)}{a^{(n1)/2}}, \quad n -1,\ a 0$$这个式子看起来吓人但它是母公式上面表格里每一行都是它的特例。比如 $n0$$\frac{1}{2}\Gamma(1/2)/a^{1/2} \frac{1}{2}\sqrt{\pi}/a^{1/2}$对上了$n2$$\frac{1}{2}\Gamma(3/2)/a^{3/2} \frac{1}{2}\cdot\frac{\sqrt{\pi}}{2}/a^{3/2}$也对上了。我建议的做法是只背母公式再记住 $\Gamma(1/2)\sqrt{\pi}$ 和 $\Gamma(n1)n\Gamma(n)$ 这两条递推关系剩下的现场推三十秒的事比记七行表格可靠得多。2. 一台引擎带动三个学科物理、统计、数学里的同一种结构2.1 大学物理麦克斯韦速率分布是怎么配出来的麦克斯韦速率分布的密度函数是$$f(v) 4\pi\left(\frac{m}{2\pi k T}\right)^{3/2} v^2 \exp!\left(-\frac{mv^2}{2kT}\right)$$这个式子里有两个东西需要理解。前面那个 $4\pi v^2$ 是速度空间里半径为 $v$ 的球面面积因为速度是三维矢量所有大小相同、方向不同的速度都落在这个球面上所以速率为 $v$的概率正比于这个球面积。后面的指数项来自能量分布 $e^{-E/kT}$其中 $E \frac{1}{2}mv^2$。归一化常数 $4\pi(m/2\pi kT)^{3/2}$ 就是为了让 $\int_0^\infty f(v),dv 1$ 恰好成立。验证一下这个归一化令 $a m/(2kT)$则积分是 $4\pi(m/2\pi kT)^{3/2}\int_0^\infty v^2 e^{-av^2}dv$。内层用表格里的第三行得 $\frac{\sqrt{\pi}}{4a^{3/2}}$。代回去$$4\pi\left(\frac{m}{2\pi kT}\right)^{3/2}\cdot\frac{\sqrt{\pi}}{4\left(\frac{m}{2kT}\right)^{3/2}} \pi\sqrt{\pi}\left(\frac{1}{\pi kT}\right)^{3/2}(kT)^{3/2} \pi\sqrt{\pi}\cdot\pi^{-3/2} 1$$正好是 1。这个验证值得亲手做一次因为它是所有归一化问题的模板把分布写成几何因子 × 指数因子 × 常数常数就是用它自己反推出来的。接下来是三个特征速率。最概然速率靠求导对 $v^2e^{-av^2}$ 求导得 $2ve^{-av^2} - 2av^3e^{-av^2}$令其为零得 $v^2 1/a$所以 $v_p \sqrt{2kT/m}$。平均速率要用 $n3$ 的那一行$\langle v \rangle 4\pi(m/2\pi kT)^{3/2}\cdot\frac{1}{2a^2}$代入 $a m/(2kT)$ 化简得到 $\sqrt{8kT/\pi m}$。方均根速率用 $n4$ 那一行得到 $\langle v^2\rangle 3kT/m$开方即 $\sqrt{3kT/m}$。三个速率之比是固定的$v_p : \langle v\rangle : v_{rms} 1 : 2/\sqrt{\pi} : \sqrt{3/2} \approx 1 : 1.128 : 1.225$。这个比例跟气体种类、温度都无关是我在读文献时判断某个作者有没有抄错的一个快速验算手段。具体数字也值得算一遍。取氮气$M 28\ \mathrm{g/mol}$$T 300\ \mathrm{K}$。用摩尔形式更省事$v_p \sqrt{2RT/M}$$R 8.314\ \mathrm{J/(mol\cdot K)}$$2 \times 8.314 \times 300 / 0.028 178157$开方得 $422.1\ \mathrm{m/s}$。平均速率 $ 422.1 \times 1.128 476.2\ \mathrm{m/s}$方均根 $ 422.1 \times 1.2247 517.0\ \mathrm{m/s}$。这三个数在空气中声速约 343 m/s和分子平均热运动速度之间建立了一个很直观的量级感声波的传播速度比分子平均速率还慢说明碰撞传递的效率是有限的。注意这里最容易犯的错是单位。$m$ 必须是单个分子质量单位是 kg如果用摩尔质量 28就必须用 $R$ 而不是 $k$。混用的话结果会差出约 31.6 倍的平方根量级一眼就能看出不对但很多人第一反应是公式抄错了反复检查公式而不是单位。2.2 概率统计正态分布从哪来为什么长相永远一样正态分布的密度是 $\frac{1}{\sqrt{2\pi}\sigma}e^{-(x-\mu)^2/(2\sigma^2)}$前面那个 $1/(\sqrt{2\pi}\sigma)$ 完全来自高斯积分。验证方式令 $u x - \mu$$\int_{-\infty}^{\infty}e^{-u^2/(2\sigma^2)}du$这里 $a 1/(2\sigma^2)$结果是 $\sqrt{\pi/a} \sqrt{2\pi\sigma^2} \sqrt{2\pi}\sigma$除以它自身正好是 1。这一步解释了为什么正态分布的归一化常数里必须有 $\sqrt{2\pi}$它不是一个约定的写法而是积分算出来的必然结果。顺着这条路往前走会碰到两个更有价值的结论。第一个是特征函数。$\varphi(t) E[e^{itX}]$对标准正态$$\varphi(t) \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty} e^{itx}e^{-x^2/2}dx$$这里用到带线性项的那个公式$a 1/2$$b it$$b^2/(4a) -t^2/2$于是结果是 $\frac{1}{\sqrt{2\pi}}\sqrt{2\pi}e^{-t^2/2} e^{-t^2/2}$。这个结果有多重要它说明正态分布在傅里叶变换下形状不变还是高斯的。所有中心极限定理的证明、所有关于独立正态变量相加仍是正态的结论根子都在这里。如果只把高斯积分当成一个孤立的反常积分去学就会错失这条线索。第二个是多维情形。$n$ 维高斯密度的归一化常数是 $1/\sqrt{(2\pi)^n\det\Sigma}$其中 $\Sigma$ 是协方差矩阵。这个 $\det$ 从哪来来自多维高斯积分公式 $\int e^{-\frac{1}{2}\mathbf{x}^T A\mathbf{x}}d^n\mathbf{x} (2\pi)^{n/2}/\sqrt{\det A}$。直观理解是协方差矩阵的行列式衡量了这个高斯椭球的体积行列式越大分布摊得越开峰值就越低所以前面的系数就越小。做机器学习的人天天在算马氏距离和高斯似然如果不知道这个系数的来源调参时会完全没有手感。2.3 高等数学Γ 函数、误差函数与反常积分的训练场高斯积分在数学分析里扮演的角色是入口。从它出发可以自然引出 Gamma 函数$\Gamma(z) \int_0^\infty t^{z-1}e^{-t}dt$令 $t x^2$ 就得到 $\Gamma(1/2) 2\int_0^\infty e^{-x^2}dx \sqrt{\pi}$。这是 Gamma 函数在非整数点上少数几个能精确求值的位置之一而且它的值直接来自高斯积分。另一个入口是误差函数及其补函数。$\operatorname{erf}(x) \frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}dt$性质是奇函数、$\operatorname{erf}(\infty)1$、$\operatorname{erf}(0)0$。补函数 $\operatorname{erfc}(x) 1 - \operatorname{erf}(x)$ 在数值上更常用原因马上会在第四节讲到当 $x$ 较大的时候直接算 $1 - \operatorname{erf}(x)$ 会因为两个相近数相减而损失精度。再往上一层是渐近展开。$\int_x^\infty e^{-t^2}dt$ 在 $x$ 很大时的主项是 $\frac{e^{-x^2}}{2x}$这个式子在做尾部概率估算时极其有用。比如想估算 $x3$ 以外的概率$\frac{e^{-9}}{6} \approx 2.06\times 10^{-5}$乘 2 得约 $4.1\times 10^{-5}$和精确值 $2.7\times 10^{-3}$这是标准正态尾部的对应值注意两者的被积函数差一个 $1/2$ 因子需要换算清楚。这个渐近公式的价值在于不用查表、不用跑代码心算就能估出数量级做建模时判断这个尾部要不要认真处理非常快。3. 手推三遍不算多配方法、矩公式与多维情形3.1 加了一个一次项怎么处理配方法实际遇到的积分经常带着线性项比如 $\int e^{-ax^2bx}dx$。很多人第一反应是分部积分那会越做越乱。正确做法是配方法跟解一元二次方程的思路一模一样$$-ax^2 bx -a\left(x^2 - \frac{b}{a}x\right) -a\left(x - \frac{b}{2a}\right)^2 \frac{b^2}{4a}$$于是积分变成 $e^{b^2/(4a)}\int e^{-a(x - b/(2a))^2}dx$。后面这个积分只是把标准高斯积分平移了一下平移不改变整个实数轴上的积分值所以等于 $\sqrt{\pi/a}$。最终结果 $\sqrt{\pi/a},e^{b^2/(4a)}$。这个式子的意义远不止一个积分技巧。特征函数 $e^{-t^2/2}$、生成函数、路径积分里的高斯型积分、乃至热传导方程的基本解全靠这一步。我在教别人时通常会说看到高斯型里出现一次项条件反射地配方不要想别的。提示配方之后那个常数因子 $e^{b^2/(4a)}$ 只能提到积分号外面绝对不能和前面的归一化系数混在一起约掉这是手算时最常见的丢项位置。3.2 矩公式把 Γ 函数当查表用$\int_0^\infty x^n e^{-ax^2}dx$ 这一族积分的推导有一个统一的换元令 $t ax^2$则 $x \sqrt{t/a}$$dx \frac{1}{2\sqrt{at}}dt$。代入$$\int_0^\infty \left(\frac{t}{a}\right)^{n/2}e^{-t}\cdot\frac{1}{2\sqrt{at}}dt \frac{1}{2a^{(n1)/2}}\int_0^\infty t^{(n-1)/2}e^{-t}dt \frac{\Gamma!\left(\frac{n1}{2}\right)}{2a^{(n1)/2}}$$这个推导只有三行但把整族公式一次性解决了。我建议每个人都自己推一遍推完之后这些公式就再也不会记错了。之后需要用具体的 $n$只需查 Gamma 函数的两个性质$\Gamma(1/2) \sqrt{\pi}$$\Gamma(n1) n\Gamma(n)$。举两个实操中一定会用到的例子。求 $\int_0^\infty x^2 e^{-ax^2}dx$$n2$$\Gamma(3/2) \frac{1}{2}\Gamma(1/2) \frac{\sqrt{\pi}}{2}$结果 $\frac{\sqrt{\pi}}{4a^{3/2}}$。求 $\int_0^\infty x^4 e^{-ax^2}dx$$n4$$\Gamma(5/2) \frac{3}{2}\cdot\frac{1}{2}\sqrt{\pi} \frac{3\sqrt{\pi}}{4}$结果 $\frac{3\sqrt{\pi}}{8a^{5/2}}$。这两个公式分别对应二阶矩和四阶矩也就是方差和峰度计算中反复出现的东西。还有一个容易忽略的边界条件母公式要求 $n -1$ 且 $a 0$。$n -1$ 时积分是对数发散的$a \le 0$ 时积分在无穷远处发散。写代码做数值积分时如果参数取到了这些边界quad会给出 warning 或者返回 inf这时候先检查参数范围别急着怪库函数。3.3 多维高斯行列式是怎么冒出来的多维情形是很多人真正需要但最不敢碰的部分。设 $A$ 是对称正定矩阵考虑 $J \int_{\mathbb{R}^n} e^{-\frac{1}{2}\mathbf{x}^TA\mathbf{x}}d^n\mathbf{x}$。做法是正交对角化因为 $A$ 对称存在正交矩阵 $Q$ 使 $Q^TAQ \Lambda \mathrm{diag}(\lambda_1,\dots,\lambda_n)$。令 $\mathbf{x} Q\mathbf{y}$由于 $Q$ 正交雅可比行列式绝对值为 1于是$$J \int e^{-\frac{1}{2}\sum_i \lambda_i y_i^2}d^ny \prod_{i1}^{n}\int_{-\infty}^{\infty}e^{-\lambda_i y_i^2/2}dy_i \prod_i\sqrt{\frac{2\pi}{\lambda_i}} \frac{(2\pi)^{n/2}}{\sqrt{\det A}}$$关键在最后一步$\prod_i \lambda_i \det\Lambda \det A$行列式的乘积性质在这里直接变成了开方后的分母。这个结果解释了两件事第一为什么多维高斯的归一化常数是 $1/\sqrt{(2\pi)^n\det\Sigma}$第二为什么 $\det\Sigma$ 必须为正——协方差矩阵半正定如果不满秩比如两个特征完全共线行列式为零归一化常数发散这说明该分布在 $n$ 维空间里塌到了一个低维子空间上不存在关于 $n$ 维勒贝格测度的密度。实际代码里踩过的一个坑协方差矩阵是从数据估计出来的浮点误差可能让它出现微小的负特征值$\det$ 变成负数开方得到 NaN。处理办法是加一个小小的对角正则项 $\Sigma \epsilon I$或者在求 $\det$ 前先做对称化和特征值截断。这个细节在标准教材里基本不会提但做工程实现时几乎必然会撞上。4. 代码里怎么落地数值积分、蒙特卡洛与精度控制4.1 先用 quad 和 gamma 互相对账写数值代码的第一个建议是不要让公式和代码同时出错。做法是用两条完全不同的路径算同一个量然后对比。一条路径是数值积分另一条是解析公式。下面这段我常用的对账脚本import numpy as np from scipy.integrate import quad from scipy.special import gamma a 1.7 for n in [0, 1, 2, 3, 4]: num, err quad(lambda x: x**n * np.exp(-a * x**2), 0, np.inf) ana 0.5 * gamma((n 1) / 2) / a**((n 1) / 2) print(fn{n} num{num:.12f} ana{ana:.12f} rel{abs(num-ana)/ana:.2e})我这边跑出来的结果是每一行的相对误差都在 $10^{-15}$ 量级或更小quad返回的 error 估计也在 $10^{-10}$ 以下。这就说明公式和代码至少有一个是对的通常这就够了。如果某一行对不上先单独检查那一行的解析式因为quad对这类光滑、快速衰减的被积函数几乎不会出错。值得强调的一点quad处理无穷上限时会自动做变量替换把 $[0,\infty)$ 映射到有限区间所以直接传np.inf是正确用法不需要手动截断。反过来如果手动截断到 $[0, 40]$对 $a1$ 的情况误差可以忽略但对 $a 10^{-4}$ 的情况积分尺度在 $1/\sqrt{a} 100$截断到 40 就只算了一半多一点结果会明显偏小。这个坑我见过不止一次本质上是忘了检查被积函数的有效宽度是多少。4.2 蒙特卡洛提议分布选对了方差可以接近零数值积分之外蒙特卡洛是另一条常用路径尤其是在高维情形下。标准做法是把积分写成期望$\int f(x)dx \int \frac{f(x)}{p(x)}p(x)dx E_p\left[\frac{f(x)}{p(x)}\right]$其中 $p$ 是任选的概率密度。这个式子成立的前提是 $p$ 在 $f \ne 0$ 的地方都为正。拿 $I \int_{-\infty}^{\infty}e^{-x^2}dx$ 做实验用三种提议分布对比提议分布 $p(x)$权重 $w(x) e^{-x^2}/p(x)$单样本方差一百万个样本的相对误差均匀分布 $U(-3,3)$$6e^{-x^2}$较大且有截断偏差有系统偏差约 $10^{-5}$标准正态 $N(0,1)$$\sqrt{2\pi}e^{-x^2/2}$约 0.486约 $4\times 10^{-4}$$N(0, 1/2)$恒等于 $\sqrt{\pi}$恰好为 0机器精度第三种情况值得单独说。因为 $e^{-x^2} \sqrt{\pi}\cdot\frac{1}{\sqrt{\pi}}e^{-x^2}$而 $\frac{1}{\sqrt{\pi}}e^{-x^2}$ 恰好就是 $N(0, 1/2)$ 的密度函数所以权重 $w$ 是一个常数 $\sqrt{\pi}$方差为零无论采多少样本估计值都严格等于 $\sqrt{\pi}$这里指的是方差为零这个理论事实实际中浮点误差还是有的。这说明了重要性采样的核心思想当提议分布的形状和被积函数完全匹配时估计是零方差的。这个例子看起来像是巧合其实不是。它揭示了一个通用的设计原则如果你能找到一个和被积函数形状接近、又容易采样的分布作为提议分布方差就能大幅降低。实际工作中常见的手段是先做拉普拉斯近似用一个高斯去逼近目标函数再拿这个高斯当提议分布。这在高维贝叶斯计算里是标准操作。import numpy as np rng np.random.default_rng(20240607) n 1_000_000 # 提议分布 N(0, 1/2)权重恒为 sqrt(pi) x rng.normal(0.0, np.sqrt(0.5), n) w np.exp(-x**2) / (np.exp(-x**2 / (2 * 0.5)) / np.sqrt(2 * np.pi * 0.5)) print(w.mean(), w.std()) # 1.7724538509055159 0.0我实测这段的输出就是 $\sqrt{\pi}$ 和零标准差显示为 0 是因为浮点层面所有 $w$ 完全一样。第一次看到这个结果会觉得有点作弊但它确实是用蒙特卡洛算出了高斯积分而且不需要事先知道答案——只要知道提议分布和被积函数成比例这个结构就行。4.3 截断区间、下溢与相减抵消第三节讲的是算法这里讲三个实际调试中最常见的精度问题。第一个是截断。很多人图省事把 $(-\infty,\infty)$ 截成 $[-5,5]$ 或者 $[-6,6]$。$[-5,5]$ 的截断误差是 $2\int_5^\infty e^{-x^2}dx$用渐近公式估一下$\frac{e^{-25}}{10}\approx 1.4\times 10^{-12}$对双精度来说已经足够好。但如果截到 $[-3,3]$误差是 $2\cdot\frac{e^{-9}}{6}\approx 4.1\times 10^{-5}$相对误差 $2.3\times 10^{-5}$只有四五位有效数字。问题在于这个误差是系统性的不会随着采样点变多而消失很多人把它误当成随机误差调了半天样本量也没用。我的一般做法是对 $e^{-x^2}$ 这类直接写np.inf必须截断时先按 $\sqrt{-\ln \epsilon}$ 估一下该截多宽$\epsilon$ 取 $10^{-16}$ 时约等于 6所以 $[-6,6]$ 是个安全的经验值。第二个是下溢。np.exp(-x**2)在 $x 30$ 时是 $e^{-900}\approx 10^{-391}$双精度直接下溢成 0。在对数似然、对数配分函数的计算中如果先把很多这样的因子乘起来再取对数会得到 $-\infty$ 或者 0整个结果就废了。标准做法是全程在对数域操作用logsumexp之类的技巧。这一点在统计建模里几乎是必修课。第三个是相减抵消最典型的是尾部概率。求 $P(Z 8)$直观写法是1 - norm.cdf(8)。但norm.cdf(8)约等于 $0.9999999999999993$用 1 去减只剩下 $6.7\times 10^{-16}$ 这个量级的差值而真实值是 $6.22\times 10^{-16}$相对误差接近 7%。正确写法是用生存函数或者对数尾概率import numpy as np from scipy.stats import norm from scipy.special import log_ndtr print(1 - norm.cdf(8)) # 6.66e-16 有效位损失严重 print(norm.sf(8)) # 6.22e-16 正确 print(np.exp(log_ndtr(-8))) # 6.22e-16 对数域极端尾部也稳log_ndtr的好处是它在 $x$ 特别大时也稳定做多重检验校正或者极端事件概率估算时会用到。这条经验适用于所有用 1 减去一个接近 1 的数的场景凡是看到这种结构就应该换一个函数。5. 常见问题速查六个高频错误与排查路径5.1 错误对照表下面这张表是我这几年在代码审查和答疑里反复见到的按出现频率排序症状典型原因排查办法归一化常数差 $\sqrt{2}$把 $\int e^{-x^2/2}dx$ 当成 $\sqrt{\pi}$认清指数里的系数$e^{-ax^2}$ 对应 $\sqrt{\pi/a}$结果差 $2\pi$ 或 $\pi$极坐标换元漏了雅可比因子 $r$回看二重积分那一步检查 $dx,dy \to r,dr,d\theta$速率算出来差 30 多倍摩尔质量配了 $k$或者分子质量配了 $R$统一用 $RT/M$ 或 $kT/m$不要混搭数值积分结果偏小且不随样本增加改善积分区间截断过窄按 $\sqrt{-\ln\epsilon}$ 估算有效宽度协方差矩阵开方后是 NaN$\det\Sigma$ 因浮点误差为负对称化 加 $\epsilon I$ 正则尾部概率偏差大$1 - \operatorname{cdf}$ 相减抵消改用sf或log_ndtr这张表里的前两条占了绝大多数。尤其是第一条$e^{-x^2}$ 和 $e^{-x^2/2}$ 长得太像一忙起来就混。我的习惯是每次写归一化之前先在纸上把指数里的系数明确标成 $a$然后把 $a$ 的表达式写在旁边最后再代回去。这个动作只花十秒钟但能省掉半小时的调试。5.2 几条个人避坑习惯第一条习惯是先量纲再数字。任何高斯型的公式先检查结果的数量级和量纲对不对。比如速率分布的归一化结果必须是无量纲的 1如果是 $\mathrm{s/m}$那一定哪里错了。这个检查不需要算具体数值扫一眼单位就能发现大部分错误。第二条是至少两条路径对账。公式、数值积分、蒙特卡洛、查表任选两条互相验证。上文那个quad和gamma对账的脚本我基本上是复制粘贴到每个新项目里的五十行以内的成本换的是心安。第三条是参数范围显式断言。高斯积分里 $a0$ 是个硬条件代码里加一行assert a 0的成本几乎为零但能拦住从数据里估计出来的负方差、零特征值等问题。我在一个拟合项目里就是因为没加这个断言$\det$ 为负的样本悄悄产生了 NaN最后传播到整个损失函数花了一整天才定位到。第四条是不要在线性化之后做高斯拟合。这是拟合高斯曲线时的经典陷阱为了把非线性最小二乘变成线性的对数据取对数然后对 $\ln y$ 做二次拟合。问题在于取对数会改变噪声的分布小 $y$ 处的相对误差被放大而高斯尾部恰恰是 $y$ 最小的地方结果拟合出来的 $\sigma$ 系统性偏小。正确的做法是直接在原空间做非线性最小二乘或者用最大似然。这个坑和积分本身没直接关系但凡是处理高斯型数据的人都会遇到值得单独记一笔。6. 最后聊几句使用习惯我现在处理高斯积分的习惯跟刚学的时候完全反过来了。以前是遇到一个积分就去翻公式表找一个长得最像的然后照着套现在是先看指数里是什么把系数提出来命名为 $a$把线性项命名为 $b$然后直接写母公式剩下的交给 Gamma 函数。这个转变带来的最大好处不是算得快而是不容易错——因为母公式只有一个需要记住的边界条件也只有 $a0$、$n-1$ 这两条检查面小了很多。另一个体会是数值和解析一定要配对使用。纯解析容易在系数上翻车纯数值容易在截断和精度上翻车两者互相印证才能既快又稳。我现在的默认流程是先用解析式算出闭式结果再用quad在无穷区间上验一遍误差超过 $10^{-10}$ 就停下来查如果这个积分要在