高斯积分这个名字，很多人是在大学物理课上第一次听到的，但真正被它咬一口，往往是在写代码或者对公式的时候。前阵子一个做分子运动模拟的同学拿着速率分布曲线来找我，说怎么画都归一化不到 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^2+y^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^2+y^2)}$ 这个函数具有旋转对称性，在直角坐标下两个变量纠缠不清，在极坐标下直接拆成了两个独立的部分，内层关于 $r$ 的积分因为多了那个 $r$ 因子，正好可以用换元 $u = r^2$ 一步做出来。如果少了极坐标里那个 $r$，看起来形式一样，但内层积分就变成了 $\int e^{-r^2}dr$，又回到原点。
注意：这个技巧只对"整个实数轴"成立。如果换成 $\int_0^\infty$，用对称性可以直接得到半个结果 $\sqrt{\pi}/2$；但如果上限是有限数，比如 $a$，那就没有闭式解了，只能用误差函数表示。
我个人的体会是，带学生推这一步时，最值钱的地方不是结果，而是让他们自己意识到"什么时候该考虑极坐标"。判据很简单：被积函数里出现了 $x^2+y^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^2+bx}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{n+1}{2}\right)}{a^{(n+1)/2}}, \quad n > -1,\ a > 0$$
这个式子看起来吓人，但它是"母公式"，上面表格里每一行都是它的特例。比如 $n=0$：$\frac{1}{2}\Gamma(1/2)/a^{1/2} = \frac{1}{2}\sqrt{\pi}/a^{1/2}$，对上了；$n=2$：$\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(n+1)=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}$。平均速率要用 $n=3$ 的那一行：$\langle v \rangle = 4\pi(m/2\pi kT)^{3/2}\cdot\frac{1}{2a^2}$，代入 $a = m/(2kT)$ 化简得到 $\sqrt{8kT/\pi m}$。方均根速率用 $n=4$ 那一行，得到 $\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}$，这个式子在做尾部概率估算时极其有用。比如想估算 $x=3$ 以外的概率，$\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^2+bx}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^{(n+1)/2}}\int_0^\infty t^{(n-1)/2}e^{-t}dt = \frac{\Gamma!\left(\frac{n+1}{2}\right)}{2a^{(n+1)/2}}$$
这个推导只有三行，但把整族公式一次性解决了。我建议每个人都自己推一遍，推完之后这些公式就再也不会记错了。之后需要用具体的 $n$，只需查 Gamma 函数的两个性质：$\Gamma(1/2) = \sqrt{\pi}$，$\Gamma(n+1) = n\Gamma(n)$。
举两个实操中一定会用到的例子。求 $\int_0^\infty x^2 e^{-ax^2}dx$：$n=2$，$\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$：$n=4$，$\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_{i=1}^{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(f"n={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]$，对 $a=1$ 的情况误差可以忽略，但对 $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对账的脚本我基本上是复制粘贴到每个新项目里的，五十行以内的成本，换的是心安。
第三条是"参数范围显式断言"。高斯积分里 $a>0$ 是个硬条件，代码里加一行assert a > 0的成本几乎为零，但能拦住从数据里估计出来的负方差、零特征值等问题。我在一个拟合项目里就是因为没加这个断言，$\det$ 为负的样本悄悄产生了 NaN，最后传播到整个损失函数，花了一整天才定位到。
第四条是"不要在线性化之后做高斯拟合"。这是拟合高斯曲线时的经典陷阱：为了把非线性最小二乘变成线性的，对数据取对数，然后对 $\ln y$ 做二次拟合。问题在于取对数会改变噪声的分布，小 $y$ 处的相对误差被放大，而高斯尾部恰恰是 $y$ 最小的地方，结果拟合出来的 $\sigma$ 系统性偏小。正确的做法是直接在原空间做非线性最小二乘，或者用最大似然。这个坑和积分本身没直接关系，但凡是处理高斯型数据的人都会遇到，值得单独记一笔。
6. 最后聊几句使用习惯
我现在处理高斯积分的习惯，跟刚学的时候完全反过来了。以前是遇到一个积分就去翻公式表，找一个长得最像的，然后照着套；现在是先看指数里是什么，把系数提出来命名为 $a$，把线性项命名为 $b$，然后直接写母公式，剩下的交给 Gamma 函数。这个转变带来的最大好处不是算得快，而是不容易错——因为母公式只有一个，需要记住的边界条件也只有 $a>0$、$n>-1$ 这两条，检查面小了很多。
另一个体会是，数值和解析一定要配对使用。纯解析容易在系数上翻车，纯数值容易在截断和精度上翻车，两者互相印证才能既快又稳。我现在的默认流程是：先用解析式算出闭式结果，再用quad在无穷区间上验一遍，误差超过 $10^{-10}$ 就停下来查；如果这个积分要在
