Monte Carlo · Buffon · PRNG

随机数与随机模拟：蒙特卡洛方法

概率难算的时候还有一条退路：把试验重复很多次，用频率去顶替概率。这条退路的合法性来自大数定律，见频率与概率：用数据估计可能性；而它的代价写在误差里，误差按 $1/\sqrt n$ 收缩，慢到成为整套方法的主要约束。把这条思路用在几何概型上就得到蒙特卡洛方法（Monte Carlo method）：先把待求的量写成某个区域的测度之比，再撒点数比例。历史上第一个这样的实验是 1777 年的布丰投针（Buffon's needle），它算出的常数里带着 $\pi$ 。所有这些都建立在计算机能产出「随机」数之上，而 Math.random() 产出的并不随机，只是难以预测。

1 · 频率估计的误差速率

SE = √(p(1−p)/n)

把一次试验重复 $n$ 次，事件发生 $k$ 次，用 $\hat p_n = k/n$ 估计 $p$ 。 $k$ 服从二项分布 $B(n, p)$ ，于是 $\hat p_n$ 的期望恰是 $p$ （无偏），方差是 $p(1-p)/n$ 。它的平方根

\mathrm{SE}(\hat p_n) = \sqrt{\frac{p(1-p)}{n}}

称为标准误（standard error），是估计值典型偏离真值的幅度。分子与 $n$ 无关，分母是 $\sqrt n$ ，所以精度与样本量之间是一条很陡的兑换比：样本量翻 4 倍，误差减半；想让小数点后多出一位有效数字，样本量要乘 100。这一条决定了随机模拟能做什么、不能做什么。

注 · 「误差按 $1/\sqrt n$ 收缩」是一句关于典型幅度的陈述，不是关于单次结果的保证。任何一次模拟都可能撞上大偏差，只是概率随 $n$ 迅速变小。中心极限定理给出更细的口径： $n$ 较大时 $\hat p_n$ 近似服从以 $p$ 为均值、 $\mathrm{SE}$ 为标准差的正态分布，落在 $p \pm 2\,\mathrm{SE}$ 内的概率约 95%。

2 · 蒙特卡洛方法

π̂ = 4 × 命中数 / 撒点数

在边长 2 的正方形里随机撒点，落进其内切单位圆的概率是面积比 $\pi/4$ 。反过来，用命中比例乘 4 就得到 $\pi$ 的一个估计 $\hat\pi = 4k/n$ 。它的标准差是 $4\,\mathrm{SE}$ ，把 $p = \pi/4$ 代进去得 $4\sqrt{p(1-p)} \approx 1.6422$ 除以 $\sqrt n$ ： $n = 10^6$ 时约 0.00164，即百万个点只买到小数点后两位半。

用 core/probstat.ts 的 montePi 各撒 100 万点，seed 取 1、2、3 时 $\hat\pi$ 分别是 3.140588、3.143188、3.142116，与 $\pi$ 相差 0.001005、0.001595、0.000523，都在预测的 0.001642 上下。把每个 $n$ 重复 1000 次算均方根误差再乘 $\sqrt n$ ，理论值应恒为 1.6422，实测（mulberry32 / Math.random()）在 $n = 10^3$ 到 $10^6$ 之间是 1.6497 / 1.6180、1.6589 / 1.5437、1.6261 / 1.6725、1.7089 / 1.6910。四个数量级上没有系统性偏离，速率就是 $1/\sqrt n$ ，一点没多也一点没少。

同一套做法可以估任何能判定「点在不在里面」的区域面积：把区域装进一个已知面积的矩形，命中比例乘矩形面积即得。写成积分的形式是 $\int_a^b f \approx (b-a)\cdot \frac{1}{n}\sum_i f(U_i)$ ，其中 $U_i$ 在 $[a, b]$ 上均匀。用 monteIntegrate 取 100 万个样本算 $\int_0^1 \sqrt{1-x^2}\,\mathrm{d}x$ 得 0.785743，真值 $\pi/4 = 0.785398$ 。

图 2-1 · 左侧是撒点估

\pi

的过程，右侧是误差随撒点数的双对数曲线与斜率

-1/2

的参考线。可持续撒点，观察误差曲线贴着参考线下行。

一维、二维的定积分不该用蒙特卡洛：梯形法误差是 $O(h^2)$ 、辛普森法 $O(h^4)$ ，同样的求值次数下精度高出好几个数量级。蒙特卡洛的价值在高维——它的 $1/\sqrt n$ 速率与维数无关，而网格法要在 $d$ 维上保持同样的网格密度，点数随 $d$ 指数增长。二十维的积分没有别的通用办法。

3 · 布丰投针

P = 2L / (πd)

地板上画一族间距 $d$ 的平行线，把一根长 $L \le d$ 的针随机扔上去，问它与某条线相交的概率。

一次投掷的结果由两个量确定：针心到最近一条线的距离 $y$ ，以及针与线族的夹角 $\theta$ 。 $y$ 在 $[0, d/2]$ 上均匀， $\theta$ 在 $[0, \pi/2]$ 上均匀，两者独立，样本空间就是一个面积 $\pi d/4$ 的矩形，又一次落回二维面积模型，见几何概型：以测度之比定概率。针的一端到最近那条线的投影长度是 $(L/2)\sin\theta$ ，相交当且仅当 $y \le (L/2)\sin\theta$ 。这块区域的面积是

\int_0^{\pi/2} \frac{L}{2}\sin\theta \,\mathrm{d}\theta = \frac{L}{2}, \qquad P = \frac{L/2}{\pi d/4} = \frac{2L}{\pi d}.

概率里出现了 $\pi$ ，把式子倒过来便得到 $\hat\pi = 2Ln/(dk)$ ：投针可以量圆周率。 $L = d$ 时理论频率是 $2/\pi = 0.636620$ ；用 buffonFrequency 投 100 万次，三个 seed 给出 0.636089、0.636702、0.636473，反解出的 $\hat\pi$ 是 3.144214、3.141187、3.142317。

图 3-1 · 布丰投针的落针图与实时统计，压线的针高亮显示。可调针长与线距之比，对照相交频率与理论值

2L/(\pi d)

。

警示 · 投针估 $\pi$ 的效率比撒点更差： $\pi$ 出现在分母，频率的相对误差被放大，且 $L = d$ 时每次投掷只贡献不到一比特的信息。历史上 Lazzarini 在 1901 年报告用 3408 次投针得到 $355/113 = 3.1415929$ ，六位小数全对。这个精度按标准误算需要上亿次投掷，实际上是他选定了 $L/d$ 的比值并挑了一个恰好凑出答案的停止时机。投针的价值在于它是第一个把几何概率与一个超越常数连起来的实验，不在于它是求 $\pi$ 的好办法。

4 · 伪随机数与可复现性

xorshift128+ · lcg · mulberry32

Math.random() 返回的不是随机数，而是一个确定性递推产生的伪随机数（pseudorandom number）：内部有一份状态，每次调用把状态推进一步再折算成 $[0, 1)$ 里的小数。V8 用的是 xorshift128+，128 位状态、周期 $2^{128}-1$ （据 V8 团队 2015 年公布的实现，核对于 2026-08）。它统计性质够好，但 JavaScript 语言层面既不能设定它的种子，也读不到它的状态。

「不可播种」这句话经实测需要加限定。V8 有一个启动开关 --random-seed，node v26.6.0 下 node --random-seed=42 -e "console.log(Math.random())" 连跑两次都给出 0.6159659788862542。所以不可播种的准确含义是：语言层面没有对应 API，浏览器里也没有这个开关，而进程级的开关一次只能定一条流，同一页里并存多个各自可复现的实验仍然做不到。本系列的所有模拟因而自带发生器，正文引用的每个数字都能按 seed 复现。

core/probstat.ts 里有两个。lcg 是线性同余发生器，状态一个 32 位整数，递推 $s \leftarrow (1664525\,s + 1013904223) \bmod 2^{32}$ ；它是 $2^{32}$ 上的一个置换，走完整个周期才回到起点。mulberry32 是计数器式的：状态每次加常数 0x6D2B79F5，再经三轮移位与乘法混合。

两者的缺陷各不相同，而且都是跑出来才看清的。模 $2^k$ 的线性同余发生器有一条老毛病：第 $k$ 低位的周期只有 $2^k$ 。lowBits(lcg(12345), 24, 1) 打印出来是 010101010101010101010101，最低位严格交替；取低两位则是 0321 四个一循环。拿这一位去抛硬币，得到的将是正反正反的完美交替。高位没有这个问题，s / 2^32 用上了全部 32 位，20 万次抽取的均值贴在 0.5 上。mulberry32 的低位没有可见规律，代价在另一处：原以为「计数器加混合」是个双射、一个周期内不会重复取值，实测 200 万次抽取撞出 1009 个重复值（seed 换成 2、3、2026、20260818 时是 1037、1070、1056、1060），而同样 200 万次的 lcg 一个重复都没有。若把它的输出映射当成随机函数，生日问题给出的预测是 466 个，实测约 2.2 倍，说明这个混合函数的像分布并不平坦。重复率在百万分之几的量级，对估概率无碍，但那个「必是双射」的推断是错的。

图 4-1 · 两个发生器连续 32 次输出的最低若干位。可切换取位数与 seed，观察 lcg 的低位周期恒为 2 的幂，而 mulberry32 看不出周期。

建议 · 选发生器只需分清用途。要可复现的模拟与可核对的正文数字，就自带一个可播种的发生器，并把 seed 写进代码；要抽签、抽奖这类不容预测的场合，Math.random() 与自带发生器都不合格，该用 crypto.getRandomValues()。另外，不要把任何发生器的低位单独拿出来用，需要一个二值结果时比较 $u < p$ ，而不是取 $u$ 的某一位。

5 · 参考文献

Buffon, G.-L. L. (1777). Essai d'arithmétique morale. In Histoire naturelle, générale et particulière, Supplément 4. Imprimerie Royale, 46–123.
Metropolis, N., & Ulam, S. (1949). The Monte Carlo method. Journal of the American Statistical Association, 44(247), 335–341.
Badger, L. (1994). Lazzarini's lucky approximation of π. Mathematics Magazine, 67(2), 83–91.