第二季 · 第 34 章 用一百万次偶然寻找答案标记本章已读大字阅读打印本章回到上次阅读位置《算法的故事》第二季 · 第一篇：数学家的秘密武器——乌拉姆、冯·诺依曼与蒙特卡洛方法的诞生第34章配图1946年，美国洛杉矶。一位三十多岁的数学家正在从一场严重的疾病中恢复。不久之前，他曾被送进医院，接受过危险的手术。疾病让他一度失去正常说话的能力，也让他的生活发生了突然的变化。康复期间，他开始玩一种名叫 Canfield 的单人纸牌游戏。洗牌、发牌、移动纸牌。一局结束，再来一局。有些牌局能够成功，有些却无论怎样努力都无法完成。作为一名数学家，他自然而然地开始思考：如果从一副52张纸牌中随机发牌，那么一局游戏最终成功的概率究竟是多少？他试图进行数学计算。但很快发现，这个问题比表面看起来复杂得多。一副纸牌拥有天文数字般的排列方式，而游戏进行过程中，每一次选择都可能改变后面所有步骤的结果。难道必须穷尽所有可能的牌局，才能知道成功的概率？一个念头渐渐出现。如果不去计算所有可能性，而是将牌反复洗好、摆出、进行游戏，再统计其中成功的比例呢？一百局。一千局。一百万局。如果试验足够多，成功比例是否会逐渐逼近真正的概率？对于当时的数学家而言，这并不是一个全新的概率思想。但如果把它交给新近出现的电子计算机，让机器模拟成千上万次复杂事件，事情就完全不同了。这位数学家叫作斯坦尼斯瓦夫·乌拉姆（Stanisław Ulam）。他后来回忆，正是1946年康复期间的纸牌问题，启发他开始探索一种新的计算方法。这种方法最终有了一个充满戏剧性的名字：蒙特卡洛方法（Monte Carlo Method）。它不试图逐一解决每一种可能，而是让大量随机试验替我们探索答案。一个源于纸牌概率的念头，最终进入物理学、计算机科学、金融工程、医学成像、人工智能和无数科学研究领域。然而，它也带来了一个意味深长的问题：如果我们始终没有计算出严格的精确答案，为什么仍然能够相信随机试验得到的结果？第一幕：一副纸牌，与一个几乎算不完的概率问题第34章配图Canfield 是一种使用标准扑克牌的单人纸牌游戏。图为规则布局参考，不是乌拉姆1946年实际使用的牌局照片。乌拉姆1909年出生在当时奥匈帝国的伦贝格，也就是今天乌克兰的利沃夫。第一次世界大战之后，这座城市成为波兰的利沃夫（Lwów）。二十世纪二三十年代，这里拥有一个活跃的数学研究群体。数学家经常聚集在咖啡馆中讨论难题。其中最著名的地点之一，是苏格兰咖啡馆（Scottish Café）。斯特凡·巴拿赫、雨果·施泰因豪斯等数学家，都会参与这些讨论。乌拉姆很早便进入这个数学世界。他对集合论、测度论和概率论尤其感兴趣。1933年，他获得博士学位。1935年，他受到冯·诺依曼的邀请，开始前往美国开展学术交流。第二次世界大战爆发后，他留在美国。1943年，他参与了洛斯阿拉莫斯的秘密科学研究。战争结束后，他曾前往南加州大学任职。但1946年，一场严重的脑部疾病几乎中断了他的研究生涯。根据数学史资料，乌拉姆因脑炎接受了手术，康复过程对他的生活和工作都产生了影响。[MacTutor History of Mathematics][MacTutor History of Mathematics]一道被52张纸牌放大的难题我们先不讨论 Canfield 复杂的具体规则。只考虑一个更基础的问题。把52张不同的扑克牌充分洗乱，究竟可能出现多少种排列？第一张牌有52种选择。第二张有51种。第三张有50种。一直到最后一张，只剩下一种选择。因此，排列总数是：52!=52×51×50×⋯×152!=52\times51\times50\times\cdots\times152!=52×51×50×⋯×1约等于：8.07×1067\boxed{8.07\times10^{67}}8.07×1067这是一个拥有68位十进制数字的巨大整数。当然，求某种纸牌游戏的获胜概率，并不意味着一定要枚举全部 52!52!52! 种排列。数学家可以利用对称性、递推、动态规划和其他概率技巧，大幅减少计算工作。但对于规则复杂、包含许多中途状态和可能选择的纸牌游戏，直接精确分析仍可能非常困难。而且，同一发牌布局的结果可能取决于玩家采取什么策略。如果没有固定游戏策略，甚至连“成功概率”究竟指什么，都需要先明确。乌拉姆最初面对的，正是这种组合复杂性。他后来在自传《一个数学家的历险》中回忆，自己在1946年康复期间玩单人纸牌，并思考 Canfield 游戏的成功概率。他意识到，与其试图计算所有复杂的组合情况，不如实际进行大量试验，统计成功的比例。[MacTutor History of Mathematics][OSTI]一次并不产生精确答案的计算假设我们设计一种规则明确、每次从同样分布独立发牌的纸牌游戏。让程序运行1000局。其中有143局成功。那么我们可以估计：p^=1431000=0.143\hat p=\frac{143}{1000}=0.143p^=1000143=0.143也就是14.3%。这里的143只是教学假设，不是乌拉姆的原始试验数据，也不是 Canfield 的真实胜率。我们没有得到严格的精确概率。因为另外1000局，可能只有131局成功。再换一批，也可能有155局。但是只要每次试验符合规定的概率模型，重复试验便能够提供有关真实胜率的信息。而且，随着独立样本不断增加，估计的随机波动通常会减小。这正是蒙特卡洛方法最重要的出发点。我们可以不必知道每种情况的答案，只需要能够可靠地生成许多符合概率模型的情况，并观察结果。从纸牌到科学研究乌拉姆很快意识到，纸牌问题并不是孤立的。在物理系统里，也存在大量由随机事件构成的复杂过程。例如，一个粒子在某种物质中运动。它可能继续前进。可能发生散射。可能被吸收。或者触发其他物理过程。这些事件会不断产生新的分支。如果每条分支都拥有不同的发生概率，那么整个物理系统就像一场极其复杂的纸牌游戏。区别只是：纸牌游戏中的状态是扑克牌的排列。物理模拟中的状态可能是粒子的位置、方向、能量和运动历史。这两个看起来毫不相关的问题，竟然拥有相似的数学结构。都存在大量可能路径。都受到概率规律控制。都可以通过重复模拟获得统计信息。这一联系，为后来蒙特卡洛方法的计算机实现提供了关键的思想方向。历史说明：乌拉姆的纸牌回忆属于当事人的事后自述，有明确文献来源；开篇的康复环境、发牌动作和思考过程则经过文学化组织，并非有逐时记录的现场实录。第二幕：在乌拉姆之前，人类就曾尝试用偶然寻找答案第34章配图布丰（Georges-Louis Leclerc, Comte de Buffon，1707—1788），法国博物学家与数学思想家。他提出的投针问题是几何概率中的经典案例。如果我们把故事停在1946年，很容易产生一个误解：似乎在乌拉姆玩纸牌之前，没有人想过通过随机试验解决数学问题。事实并非如此。至少在十八世纪，数学家已经开始将随机实验与几何问题联系起来。其中最著名的一个例子，来自法国的布丰。把一根针扔到地板上想象地板上画着许多条平行直线。相邻两条线之间的距离都是 ddd。你拿起一根长度为 lll 的针，其中：l≤dl\le dl≤d然后，把针随机投到地板上。现在提出一个问题：这根针与某条平行线相交的概率是多少？乍看之下，这是一个很奇怪的数学问题。针可以落在不同位置。也可以形成不同角度。如果每次的位置和方向都不一样，概率应该怎样计算？布丰研究了这个问题，并在1777年出版的作品中给出了相关分析。在针的中心位置与方向满足适当的均匀随机假设时，相交概率为：P=2lπd\boxed{P=\frac{2l}{\pi d}}P=πd2l如果针长恰好等于平行线之间的距离：l=dl=dl=d就得到：P=2πP=\frac2\piP=π2这意味着，只要反复投针、统计相交比例，就能够反过来估计圆周率。[PMC][mathematicskit]图34-1：布丰投针实验相交未相交相交查看 SVG 配图原件34-004-c2ec1db405e1.svg · 点击查看大图教学示意：橙色针跨越平行线，蓝色针没有跨越。实际实验要求针的落点及角度符合设定的概率分布。这里的数学关系很令人惊叹。圆周率原本是圆周长与直径的比值。而现在，我们竟然可以借助一根针和一组平行线，通过概率估计它。不过，乌拉姆时代的方法，与布丰投针之间仍有重要区别。布丰的实验规模相对有限，属于几何概率和随机实验的早期探索。乌拉姆及其合作者所推动的，则是将随机抽样系统地用于复杂数值计算，并利用高速电子计算机处理大量模拟过程。因此，更准确的历史说法是：蒙特卡洛方法不是从无到有地发明随机抽样，而是在电子计算时代，将已有的概率思想发展成一种通用的计算技术。费米：另一个走在时代前面的人第34章配图恩里科·费米（Enrico Fermi，1901—1954），意大利裔美国物理学家。他在中子物理研究中曾采用与后来蒙特卡洛思想相近的随机抽样方法。这段历史还有一位不可忽略的人物：恩里科·费米。费米在1930年代研究中子与物质之间的相互作用时，曾经使用过随机抽样思路，估计中子在物质中的行为。这些工作并没有在当时形成后来那种普遍传播的计算机蒙特卡洛方法体系。但尼古拉斯·梅特罗波利斯后来在自己的历史回忆中，明确提到了费米更早的独立探索。1947年，费米还设计了一种机械辅助装置，后来被称为 FERMIAC，用于帮助绘制和追踪中子运动路径。它不是现代电子计算机。但它展示了同一个思想：把复杂的粒子运动过程，转化为一次次可模拟的随机选择。[UNT Digital Library][mcnp.lanl.gov]这样看来，蒙特卡洛方法的历史并不是一条笔直的线。它有十八世纪几何概率的根源。有二十世纪早期的物理实验传统。有乌拉姆从纸牌问题中获得的启发。还需要冯·诺依曼等研究者，把数学思想转变成真正能够在电子计算机上运行的程序。接下来，故事将进入洛斯阿拉莫斯。那里，数学不再只是抽象的游戏。它开始直接参与真实的科学与工程计算。第三幕：洛斯阿拉莫斯——当物理学家面对无数条可能的道路第34章配图美国新墨西哥州洛斯阿拉莫斯的历史研究环境。第二次世界大战期间，这里成为曼哈顿计划的重要研究中心；战后继续开展核物理与数值计算研究。1946年之后，乌拉姆重返洛斯阿拉莫斯。世界大战已经结束。但许多在战争中建立起来的科学研究机构并没有消失。电子计算机、核物理、数值分析和复杂工程计算，在新的环境中继续发展。研究人员仍然面临一类极其困难的问题：一个复杂物理系统，究竟会怎样随着时间演化？考虑粒子输运问题。假设有大量中子进入某种材料。一部分中子可能被吸收。另一部分可能发生散射。还有一部分可能继续穿过材料。每一次碰撞之后，粒子的方向、能量和后续行为可能发生变化。如果我们希望知道最终有多少粒子穿过材料，就需要考虑大量可能的运动历史。直接求解为什么困难？在经典数学物理中，研究者经常使用微分方程描述系统。例如热传导、电磁场和流体运动，都可以通过相应的方程建模。但如果系统中存在许多相互作用的粒子，而且每个粒子还可能经历不同的随机事件，那么方程及其边界条件可能变得非常复杂。对于某些问题，求解完整的概率分布或相关积分方程，所需计算量甚至超过当时计算设备的能力。乌拉姆提出的思路，是把问题倒过来。与其直接求解描述整个系统的复杂数学方程，不如根据已知的物理概率规律，模拟大量可能发生的粒子运动历史。对于一条粒子轨迹：从初始状态出发。根据模型选择下一次事件。更新状态。然后继续。直到这条轨迹达到某个终止条件。随后，再模拟一条新的轨迹。如此重复。最终根据大量轨迹，估计研究者关心的统计量。图34-2：从一颗粒子到千万条可能轨迹材料内部 · 随机运动区域入射穿过另一条路径吸收每条路径由模型规定的随机事件共同决定查看 SVG 配图原件34-007-3b3812132db9.svg · 点击查看大图简化的粒子输运示意图，并非物理尺度准确的轨迹计算。实际模拟还涉及碰撞截面、材料属性、粒子能量等模型参数。如果让计算机模拟一万个粒子我们可以建立一个极简的教学模型。假设某类粒子每经过一个观察阶段，会出现三种结果：20%的概率被吸收。30%的概率发生散射并进入下一阶段。50%的概率直接穿过材料。这些概率只为解释算法而设定，不对应某种真实材料的中子输运数据。计算机每次产生一个0到1之间的随机数。如果随机数小于0.2，就判定为吸收。如果落在0.2到0.5之间，就判定为散射。如果大于等于0.5，就判定为穿过。其中，散射的粒子继续进入下一轮。最终，我们就可以统计有多少粒子穿过材料，有多少粒子被吸收。这就是一个简单的随机状态转移模拟。虽然它远不及真实粒子输运模型复杂，但数学思想已经相通。一个难以直接追踪所有可能结果的系统，可以被转化成大量单独的随机实验。乌拉姆需要的不只是一个好想法在纸上模拟十条轨迹，或许还不困难。模拟一百条，就会变得繁琐。模拟一百万条，几乎无法依靠人工完成。而且，每条轨迹还可能包含许多连续的随机事件。这时，乌拉姆的想法必须与另一项新技术结合：电子计算机。1947年，冯·诺依曼、罗伯特·里希特迈尔等研究者开始推动将随机模拟用于实际计算。冯·诺依曼在1947年写给里希特迈尔的信中，描述了如何组织计算机执行有关粒子输运的随机计算。相关技术史研究将这份信视为早期电子计算机蒙特卡洛程序设计的重要文献。[计算机科学系][mcnp.lanl.gov]伪随机数在这里找到了真正的舞台这也解释了上一章的重要性。计算机需要不断进行随机选择。但是，如果每次都依赖人工掷骰子或读取外部随机数据，模拟就很难高效运行。通过伪随机数生成器，计算机可以快速获得一串用于模拟的数字。对于指定的种子和程序，这串数字还能够重复。所以，第三十三章中看似矛盾的“确定性制造随机性”，在这里产生了实际价值。一个确定性的程序，可以通过伪随机数，模拟具有不确定性的物理过程。但研究者必须确保随机数生成器的统计质量符合模拟要求，否则数字序列中的相关性和偏差，就可能污染物理结果。这已经不是一场纸牌游戏。这是数值计算的一种新范式。第四幕：ENIAC、梅特罗波利斯，以及一个来自赌场的名字第34章配图ENIAC 大型电子计算机。它不仅用于早期弹道计算，也成为战后开展随机数值模拟的重要平台。图片为设备历史照片，并非某次蒙特卡洛运行记录。1940年代的电子计算机，与我们今天熟悉的机器差别巨大。它们没有现代图形桌面。没有可以随手打开的 Python 解释器。程序的组织、输入、输出和存储，都要受到硬件能力的严格限制。有时，为了让机器执行一个复杂任务，研究团队需要花费大量精力设计程序步骤，并认真规划中间数据的处理方式。在这样的技术背景下，蒙特卡洛方法开始展现真正的计算价值。一项新的计算思想，需要许多人共同完成乌拉姆提出重要思想之后，冯·诺依曼积极推动其计算机实现。里希特迈尔参与相关物理计算问题的研究。而在洛斯阿拉莫斯，从事数值计算与计算机研究的尼古拉斯·梅特罗波利斯，也成为方法发展过程中的重要人物。这些工作不是几个数学家在纸上讨论完毕，计算机就能自动运行。实际计算还需要程序设计、设备操作、测试和验证。早期程序员与计算人员的工作，同样是这段历史的一部分。到了1948年前后，使用电子计算机进行复杂随机模拟的研究已取得实际进展。现代蒙特卡洛方法，正是在这种数学、物理与计算工程的合作中逐渐成熟。[数字对象标识符][OSTI]为什么叫“蒙特卡洛”？第34章配图摩纳哥蒙特卡洛赌场，约20世纪上半叶的历史影像风格。这个与轮盘、赌博和概率联系在一起的地名，后来成为一类科学计算方法的名称。现在来到整段故事最富有戏剧性的部分。蒙特卡洛是摩纳哥著名的地区，以赌场闻名。在赌场里，轮盘旋转。骰子滚动。玩家下注。每一次结果都带有偶然性。但大量结果背后又受到概率规律约束。对于一项通过重复随机试验来解决数学问题的方法，这个名字显得十分贴切。那么是谁提出了它？梅特罗波利斯在1987年发表的回忆文章《蒙特卡洛方法的开端》中，明确表示自己建议采用这个名字。他还回忆，乌拉姆有位亲戚喜欢前往蒙特卡洛，甚至曾以需要去那里的理由向亲人借钱。这段带有幽默色彩的家庭轶事，成为名称由来的常见解释。不过，这属于梅特罗波利斯多年后的回忆。我们可以引用这份证词，却不应该将所有细节写成有同时代会议记录证明的现场对话。可以确认的是，“Monte Carlo”这一名称在1947年的洛斯阿拉莫斯内部资料中已经出现，并在1949年的公开论文中正式进入学术传播。[Taylor & Francis Online][UNT Digital Library]1949年：一篇七页的论文第34章配图尼古拉斯·梅特罗波利斯与斯坦尼斯瓦夫·乌拉姆1949年的论文《The Monte Carlo Method》，发表于《美国统计协会杂志》第44卷，第335—341页。1949年9月，《美国统计协会杂志》发表了一篇论文。作者是：Nicholas Metropolis。S. Ulam。论文标题非常直接：The Monte Carlo Method。《蒙特卡洛方法》。全文只有七页。但它讨论的不是某个孤立的掷骰子游戏，而是一类数学物理问题的通用计算方法。论文指出，对于某些复杂的微分方程、积分微分方程，以及组合概率问题，可以通过统计抽样获得实用的近似结果。它还讨论了一个非常现代的问题：高维空间中的计算。例如，在一个二十维区域里，如果每个坐标轴划分十份，规则网格可能需要检查：102010^{20}1020个格点。这对于当时的计算设备来说难以承受。但如果使用随机抽样，研究者可以先选取一定数量的点，根据这些点判断某些统计性质。当然，样本数量是否足够、结果精度如何，还需要单独分析。这项方法并不保证任意复杂问题都能轻松解决。但它提供了一条新的道路：当穷举所有状态不现实，而我们又能可靠地抽取有代表性的样本时，统计估计就可能成为可计算的替代方案。 [数字对象标识符][web.williams.edu]原论文中还有一个值得注意的思想。不同随机样本在某些问题里可以相对独立地进行计算。这意味着，未来如果拥有许多台计算机，就可以将试验分散到不同机器上执行，再汇总结果。几十年之后，并行计算、超级计算机、GPU 和大规模分布式模拟，都会不断利用这一优势。这一点对于蒙特卡洛方法极为重要。因为它提高精度的一种基本手段，就是产生更多样本。而样本越多，所需的计算资源通常也越多。第五幕：向正方形投下一百万个点，寻找圆周率第34章配图蒙特卡洛估计圆周率的经典可视化：向正方形均匀投点，统计其中落在四分之一圆内部的比例。此图为教学示意，并不是1940年代的原始科研实验。现在，我们终于可以真正动手研究蒙特卡洛方法了。对于初学者，最适合的实验之一，就是利用随机投点估计圆周率。不过有一点需要提前说明。这并不是乌拉姆和冯·诺依曼当年开展复杂物理模拟的原始案例，而是今天常见的一种教学演示。它能帮助我们在最简单的几何环境中理解蒙特卡洛方法的数学本质。第一个问题：圆和正方形有什么关系？画一个边长为1的正方形。它的面积是：Ssquare=1S_{\text{square}}=1Ssquare=1再在这个正方形中画出一个半径为1的四分之一圆。四分之一圆的面积是：Squarter=π4S_{\text{quarter}}=\frac{\pi}{4}Squarter=4π因此，圆面积与正方形面积的比值为：SquarterSsquare=π4\frac{S_{\text{quarter}}}{S_{\text{square}}}=\frac{\pi}{4}SsquareSquarter=4π这时，问题发生了变化。我们不再直接计算圆周率。而是向正方形里随机投点。只要这些点在正方形中均匀分布，一个点落入四分之一圆的概率，就等于两个图形的面积之比。因此：P(落入圆内)=π4\boxed{P(\text{落入圆内})=\frac{\pi}{4}}P(落入圆内)=4π问题已经变成了一个概率问题。第二个问题：怎样判断一个点在圆内？设随机点的坐标为：(x,y)(x,y)(x,y)其中：0≤x≤1,0≤y≤10\le x\le1,\qquad0\le y\le10≤x≤1,0≤y≤1我们知道，单位圆内部的点满足：x2+y2≤1x^2+y^2\le1x2+y2≤1所以，每生成一个点，只需要执行三步：生成随机坐标 xxx。生成随机坐标 yyy。判断：x2+y2≤1x^2+y^2\le1x2+y2≤1如果条件成立，就把它计入圆内点数量。否则，计入圆外。完成 NNN 次实验之后，设圆内点数为 KKK。那么圆内概率的估计值就是：p^=KN\hat p=\frac KNp^=NK根据刚才的面积关系：p^≈π4\hat p\approx\frac{\pi}{4}p^≈4π因此：π^=4KN\boxed{\hat\pi=4\frac KN}π^=4NK这就是整个圆周率模拟的核心公式。看起来几乎难以置信。我们不需要计算正多边形的周长。不需要推导无穷级数。只需要会产生随机点、比较数字，再统计比例。第三步：试验次数不断增加假设第一次实验投下100个点。其中80个落在圆内。那么：π^=4×80100=3.2\hat\pi=4\times\frac{80}{100}=3.2π^=4×10080=3.2与真正的圆周率相比，这个结果还有明显误差。但它已经接近3。如果投下1000个点，其中786个落在圆内：π^=4×7861000=3.144\hat\pi=4\times\frac{786}{1000}=3.144π^=4×1000786=3.144结果接近了许多。如果投下10000个点，其中7852个落在圆内：π^=4×785210000=3.1408\hat\pi=4\times\frac{7852}{10000}=3.1408π^=4×100007852=3.1408结果甚至看起来已经相当不错。需要强调，这三组数字是专门构造的教学示例，并非一次真实连续实验的记录。而且，不能据此认为试验次数每增加十倍，误差一定按固定比例减小。真正的随机抽样存在波动。我们甚至可能看到1000次试验得到的结果，比100次试验更偏离真实值。这并不违反蒙特卡洛方法的原理。因为它承诺的是统计意义上的收敛规律，而不是每增加一个样本，答案就必定变好。第六幕：为什么随机试验值得相信？大数定律的真正意义第34章配图雅各布·伯努利（Jacob Bernoulli，1655—1705）。他对大数定律的研究，为后来通过重复试验估计概率提供了重要理论基础。现在，终于来到了本章最重要的数学问题。如果蒙特卡洛方法没有逐一检查全部可能情况，为什么它仍然能够提供可靠的答案？我们先回到一个最简单的例子。假设有一枚理想且公平的硬币。正面概率是：p=12p=\frac12p=21如果只投掷一次，结果一定是正面或者反面。你无法仅凭这一次投掷，判断正面概率究竟是多少。如果投掷10次，可能有7次正面。如果投掷100次，可能有53次正面。如果投掷1000次，可能有497次正面。随着试验次数增加，正面出现的比例往往会逐渐接近0.5。这种现象背后的严格数学规律，就是大数定律。把每次试验变成一个变量对于刚才的随机投点实验，我们可以定义一个变量：Ii={1,第 i 个点在圆内0,第 i 个点在圆外I_i= \begin{cases} 1,&\text{第 }i\text{ 个点在圆内}\\ 0,&\text{第 }i\text{ 个点在圆外} \end{cases}Ii={1,0,第 i 个点在圆内第 i 个点在圆外这个变量称为指示变量。如果抽样满足独立、均匀等假设，那么：P(Ii=1)=p=π4P(I_i=1)=p=\frac{\pi}{4}P(Ii=1)=p=4π它的期望值是：E[Ii]=pE[I_i]=pE[Ii]=p我们把全部 NNN 次结果相加：K=∑i=1NIiK=\sum_{i=1}^{N}I_iK=∑i=1NIi于是：π^=4N∑i=1NIi\hat\pi=\frac4N\sum_{i=1}^{N}I_iπ^=N4∑i=1NIi这就是蒙特卡洛估计器。根据期望的线性性质：E[π^]=4N∑i=1NE[Ii]=4p=πE[\hat\pi] =\frac4N\sum_{i=1}^{N}E[I_i] =4p =\piE[π^]=N4∑i=1NE[Ii]=4p=π因此：E[π^]=π\boxed{E[\hat\pi]=\pi}E[π^]=π这说明，该估计器是无偏的（Unbiased）。所谓无偏，并不是指每一次实验都恰好得到正确答案。而是说，如果能够想象无限次重复整套实验，这些估计值的平均数等于真正的圆周率。单次实验仍然可能偏高，也可能偏低。为什么更多试验可以降低误差？由于 IiI_iIi 只可能取0或1，它的方差是：Var(Ii)=p(1−p)\operatorname{Var}(I_i)=p(1-p)Var(Ii)=p(1−p)在样本彼此独立的情况下：Var(π^)=16p(1−p)N\operatorname{Var}(\hat\pi) =\frac{16p(1-p)}{N}Var(π^)=N16p(1−p)所以标准差为：σπ^=4p(1−p)N\boxed{\sigma_{\hat\pi} =4\sqrt{\frac{p(1-p)}N}}σπ^=4Np(1−p)注意分母里的 NNN 位于根号之内。因此，典型的随机误差尺度随以下数量减小：O(1N)\boxed{O\left(\frac1{\sqrt N}\right)}O(N1)这才是蒙特卡洛方法最重要的精度规律之一。一个不太令人愉快的结论假设你希望把随机误差的典型尺度缩小到原来的十分之一。那么所需样本数不是增加十倍。而是大约增加：102=10010^2=100102=100倍。这是平方根规律带来的直接结果。独立抽样次数估计圆周率的理论标准差（约）1000.1641,0000.05210,0000.0164100,0000.00521,000,0000.00164使用 p=π/4p=\pi/4p=π/4 计算的理论标准差。这不是某次模拟的实际误差，也不能保证每次结果均落在该范围内。如果希望计算结果达到非常高的精度，直接使用这种随机投点方法，并不是一个好的选择。例如，计算圆周率到几十位小数，数学家拥有远远高效得多的确定性算法。那么，蒙特卡洛方法到底有什么优势？答案不在于它对所有问题都最快。而在于：对于某些结构复杂、维数很高、难以直接求积分或无法穷举状态的计算，随机抽样能够提供一条可控制误差的求解路线。这是一种用统计精度换取可计算性的策略。置信区间：不只给出答案，还要说明有多不确定我们还可以进一步回答：如果某次实验得到：π^=3.13\hat\pi=3.13π^=3.13那么这个结果究竟可靠吗？对于足够大的独立样本，在通常的近似条件下，可以使用：SE^(π^)=4p^(1−p^)N\widehat{SE}(\hat\pi) =4\sqrt{\frac{\hat p(1-\hat p)}N}SE(π^)=4Np^(1−p^)估计标准误差。再用正态近似，构造一个大约95%置信水平的区间：π^±1.96 SE^(π^)\hat\pi\pm1.96\,\widehat{SE}(\hat\pi)π^±1.96SE(π^)但需要特别说明，95%置信区间并不是说“在这一次计算结束以后，真实的圆周率有95%的概率落在某个固定区间内”。对于频率学派的解释，真实参数是固定的。95%指的是：如果不断重复完整的抽样和区间构造过程，按照方法构造出的区间，长期来看约有95%能够覆盖真实参数。这个解释或许有些绕。但它恰好体现了科学计算的重要精神：我们不仅要给出一个估计值，还要诚实地说明这个估计存在多大的不确定性。对于真实研究，若样本依赖关系复杂、概率极端或样本太少，上述简单的正态近似可能并不可靠，需要使用适当的其他方法。第七幕：互动实验室——亲手向圆内投下一万个点第34章配图随机投点实验的几何原理：蓝色区域对应四分之一圆，橙色区域位于正方形中、圆外。每个点都代表一次独立抽样。现在，轮到我们自己成为实验者。你拥有一个正方形。一个四分之一圆。以及一台可以不断生成随机坐标的计算机。你要做的事情只有一件：控制试验次数，观察计算出的圆周率如何变化。你会看到，估计值有时靠近真正的圆周率，有时又离它远一些。这种反复波动，并不是程序出错。它恰好是蒙特卡洛方法的重要特征。实验一：随机投点估计圆周率图34-3：蒙特卡洛圆周率实验室调整试验次数，使用不同种子比较结果。相同种子会重复得到相同的随机序列。投点次数1,00010020,000随机数种子种子 42种子 2026种子 999011为避免点重叠过密，图中最多显示前400个样本；下方数值使用你选择的全部样本计算。圆周率估计值3.224000本次绝对误差0.082407圆内点806圆外点194估计值的收敛轨迹横轴表示累计抽样次数，纵轴表示当时的估计值。曲线并不保证单调接近真实的圆周率。本次计算结果π^=4×8061000=3.22400\hat{\pi}=4\times\frac{806}{1000}=3.22400π^=4×1000806=3.22400基于正态近似的95%置信区间（示范算法）：3.1260 — 3.3220该区间是统计估计结果，不是每次实验的确定性误差保证。本实验使用固定种子的简易伪随机生成器，适合教学，不代表科学级随机性质量验证。↶ 恢复初始设置请重点观察三个现象第一个现象：当抽样只有100次时，结果可能相当不稳定。第二个现象：把抽样次数增加到10000次，估计往往更接近圆周率，但不保证每一次都比原来准确。第三个现象：使用不同的随机数种子，即使抽样次数相同，也可能得到不同估计值。这三个现象，正是蒙特卡洛方法的基本特征。实验二：用 Python 真正运行一百万次现在，让我们使用代码完成同样的任务。python复制代码 import math import random def estimate_pi(n: int, seed: int = 2026): rng = random.Random(seed) inside = 0 for _ in range(n): x = rng.random() y = rng.random() if x * x + y * y <= 1: inside += 1 estimate = 4 * inside / n error = abs(estimate - math.pi) return estimate, error for n in [100, 1000, 10000, 100000, 1000000]: pi_value, error = estimate_pi(n) print(f"实验次数: {n:>8}") print(f"估计圆周率: {pi_value:.8f}") print(f"绝对误差: {error:.8f}") print() 请注意，这里使用了：python复制代码rng = random.Random(seed)我们没有每次都创建一个完全不可重现的随机序列。而是使用固定种子，确保重复运行程序时，可以得到相同结果。这和第三十三章讨论的伪随机数直接联系起来。现代计算机上的蒙特卡洛方法，通常就是建立在可靠的伪随机数生成器之上。随机性负责探索不同可能，确定性负责让探索过程能够重复和验证。实验三：尝试一个没有简单初等反导函数的积分估计圆周率只是入门。我们再看一个稍微接近科学计算的问题。考虑积分：I=∫01e−x2 dxI=\int_0^1 e^{-x^2}\,dxI=∫01e−x2dx如果只允许使用初等函数，直接写出它的不定积分并不容易。但我们可以把它改写成一个期望：I=E[e−U2]I=E[e^{-U^2}]I=E[e−U2]这里：U∼Uniform(0,1)U\sim\operatorname{Uniform}(0,1)U∼Uniform(0,1)也就是在0到1之间均匀抽样。于是，蒙特卡洛估计就成为：I^=1N∑i=1Ne−Ui2\boxed{\hat I=\frac1N\sum_{i=1}^{N}e^{-U_i^2}}I^=N1i=1∑Ne−Ui2用 Python 实现：python复制代码 import math import random def monte_carlo_integral(n: int, seed: int = 2026): rng = random.Random(seed) total = 0.0 total_sq = 0.0 for _ in range(n): x = rng.random() y = math.exp(-x * x) total += y total_sq += y * y estimate = total / n # 样本方差，用于估计标准误差 variance = ( total_sq - n * estimate * estimate ) / (n - 1) standard_error = math.sqrt( max(0.0, variance) / n ) return estimate, standard_error value, se = monte_carlo_integral(100000) print("积分估计:", value) print("标准误差:", se) 这个例子更接近蒙特卡洛积分的通用形式。对于一维、光滑的函数，使用数值积分库中的确定性算法通常会比这种直接随机抽样更高效。然而，如果问题扩展到许多维度，或者被积函数依赖复杂的模拟过程，情况就可能改变。这正是蒙特卡洛方法在现代科学中不断获得应用的原因。第八幕：从一个圆，到拥有二十个维度的世界第34章配图高维概率与数值模拟概念示意图。真实的高维空间无法完整直接画在二维平面上，科学家通常采用投影、切片和统计摘要进行可视化。在纸上的正方形里投点，是一个二维问题。如果我们希望计算一个三维物体的体积，就可以在立方体里随机投点。如果物体的边界能够被数学条件描述，那么：落入物体内部的比例，便可以用来估计其体积与外部立方体体积的比值。再向前一步。四维空间呢？五维呢？十维呢？二十维呢？人类无法直接想象一个完整的二十维几何物体。但计算机并不一定需要把它画出来。只要能够生成二十个坐标：(x1,x2,…,x20)(x_1,x_2,\ldots,x_{20})(x1,x2,…,x20)并且判断这个点是否满足指定的数学条件，就可以开展随机抽样。一场维数引起的计算灾难假设我们使用传统规则网格。在每个维度上放置10个候选位置。一维需要：101010个点。二维需要：102=10010^2=100102=100个点。三维需要：103=100010^3=1000103=1000个点。如果来到二十维：102010^{20}1020个点。这就是高维计算里著名的维数灾难的一种表现。1000000000001维 网格点数量: 10102维 网格点数量: 1001003维 网格点数量: 100010005维 网格点数量: 1000001000008维 网格点数量: 10000000010000000010维 网格点数量: 10000000000100000000001维2维3维5维8维10维查看图表数据dimension网格点数量1维102维1003维10005维1000008维10000000010维10000000000这里的困境不是计算机速度不够快那么简单。即使每个格点只需要非常短的处理时间，格点数量也会随着维度呈指数级增长。而蒙特卡洛方法在某些高维积分问题里，却可能表现得更有优势。因为对于独立、同分布抽样，且被估计函数的方差有限的情况，样本平均值的标准误差形式为：SE=σfNSE=\frac{\sigma_f}{\sqrt N}SE=Nσf其中 σf\sigma_fσf 是被积函数在抽样分布下的标准差。这个公式中，没有直接出现“维数”这个变量。当然，这绝不意味着蒙特卡洛方法在任意高维问题中都不受维度影响。随着维数增加，函数的方差、有效样本比例、稀有事件和其他计算困难，都可能发生严重变化。但它至少说明，某些问题可以避免直接构造随维度指数膨胀的完整规则网格。这正是1949年梅特罗波利斯与乌拉姆论文重点讨论的数学优势之一。[web.williams.edu]蒙特卡洛方法并不是没有弱点它有三个经常被初学者忽略的困难。第一，罕见事件可能极难估计。如果某事件发生概率只有百万分之一，而我们只做一万次模拟，那么大多数情况下，可能一次都观察不到。这不能证明事件不可能发生。只是说明我们抽样太少。第二，抽样分布必须符合问题模型。如果我们希望在正方形内均匀投点，却生成了主要集中在中心区域的点，那么圆周率估计就会系统性偏离正确答案。随机数看起来复杂，不意味着抽样方案正确。第三，计算更多样本也需要成本。如果每条模拟路径都需要运行昂贵的物理模型，那么生成一百万个样本，仍然可能耗费大量计算资源。因此，后来研究者又进一步发展了重要性抽样、分层抽样、控制变量和其他方差缩减技术。它们的共同目标，是用更少的计算，得到更精确的统计估计。到了1953年，梅特罗波利斯与其他研究者又发表了著名的 Metropolis 算法，为后来广泛使用的马尔可夫链蒙特卡洛方法（MCMC）奠定重要基础。需要区分的是：1949年的《蒙特卡洛方法》介绍的是一类广泛的统计模拟思想。1953年的 Metropolis 算法则是一种具体的抽样方法。二者关系密切，却不是同一个算法。[计算机科学系][PMC]第九幕：一百万次试验之后，科学家真正得到的是什么？第34章配图洛斯阿拉莫斯 MANIAC 计算机及其研究人员的历史影像。20世纪50年代，更成熟的电子计算机使大规模随机模拟研究进一步发展。1950年代，洛斯阿拉莫斯又迎来了新一代电子计算机。其中一台名叫 MANIAC。这个名称是 Mathematical Analyzer Numerical Integrator and Automatic Computer 的缩写。1952年前后，MANIAC 开始在洛斯阿拉莫斯投入科研使用。它可以执行各种复杂的数值计算，随机模拟自然也是其中的重要任务。与早期依赖人工或机械设备执行大量计算的时代相比，电子计算机赋予研究者一项全新的能力：他们可以设计一个数学世界。给这个世界设置规则。规定各种事件发生的概率。然后让计算机一次又一次地模拟其中可能发生的事情。每次模拟，都是一条可能的历史。每条历史，可能都与其他历史不同。但是当无数条历史被放在一起时，某些规律就开始显现。从单个事件到整体规律一颗粒子下一次碰撞会发生什么？我们可能无法确定。但如果知道相关物理概率规律，就可以估计大量粒子最终呈现的统计结果。一场纸牌游戏最终能否成功？在固定策略、随机发牌的模型下，单次结果不一定能提前判断。但重复进行许多局，就可以估计长期胜率。一个复杂系统在许多随机条件下会怎样表现？我们也许无法提前列出全部结果。但可以通过抽样，了解一些统计特征。这就是蒙特卡洛方法所带来的根本转变。它并没有消除偶然性。它也没有告诉我们每一个事件的确定未来。它让人类学会了另一种解决问题的方式：如果无法掌握每一种可能，就通过有代表性的样本，研究整体的概率规律。一项从秘密研究走向公共科学的技术蒙特卡洛方法在核物理研究环境中迅速成长。这段历史与战争、核武器和冷战科学密切相关，不能被简单处理成一场没有背景的数学游戏。但同一种计算思想，也逐渐进入许多和平研究领域。在医学物理中，蒙特卡洛方法可以用来模拟辐射在人体组织等材料中的输运，为剂量计算和设备研究提供帮助。在金融工程中，它可以通过模拟许多可能的市场路径，估计某些金融产品或风险指标在模型假设下的分布。在计算机图形学中，它可以用于复杂光照积分，帮助计算漫反射、间接光照和路径追踪等效果。在统计学和人工智能中，蒙特卡洛方法及其扩展可以帮助估计复杂概率分布的性质，支持贝叶斯推断和其他数值计算任务。第34章配图计算机图形学随机路径采样可以帮助估计光线在场景中的传播与间接照明。第34章配图医学物理粒子输运模型可用于研究辐射剂量分布和成像系统。金融风险模拟通过不同的假设路径研究未来结果的分布，而不是预测唯一未来。分子与材料研究随机采样帮助探索复杂系统的状态分布与统计性质。这些应用并不是简单地“产生大量随机数字”就能够完成。真正有价值的蒙特卡洛模拟，至少需要三项条件。首先，有合理的数学或物理模型。其次，有能够按照模型规定概率产生样本的方法。最后，有正确的统计估计与误差分析。其中任何一个环节存在重大缺陷，都可能让计算结果失去可信度。因此，一百万次试验并不自动意味着科学上的可靠。一百万次错误建模的试验，不能代替一个正确的数学模型。而这，或许才是本章最应该记住的原则。尾声：偶然性，成为一种可以设计的计算方法让我们再次回到1946年。那时，乌拉姆正在从严重疾病中恢复。他面对的，是一副普通的扑克牌。一个看起来很小的问题：这局牌究竟有多大概率成功？如果试图把全部可能情况写在纸上，这个问题可能复杂得令人望而却步。但如果把问题改成重复试验，事情便出现了新的可能。这种思路并没有凭空创造概率论。它建立在早期数学家对随机事件、几何概率和统计规律的长期研究之上。但乌拉姆看到了一个新的机会。电子计算机的出现，使大规模随机试验第一次拥有了前所未有的实际可行性。冯·诺依曼帮助推动这种思想进入计算机。梅特罗波利斯参与了它的发展与传播，并留下了关于命名过程的重要回忆。洛斯阿拉莫斯的科学家、程序员和工程人员，又共同把这些思想变成真正可以运行的数值计算。1949年，一篇只有七页的论文，让更广阔的学术界认识了蒙特卡洛方法。此后，这套思想不断生长。从粒子的运动，到复杂积分。从物理统计，到金融风险。从医学研究，到计算机图形学。人类不再只依赖直接公式和穷举，寻找数学问题的答案。他们开始创造大量符合模型的随机世界。观察这些世界。统计其中发生的事件。再利用概率论，将有限试验的结果转化为对整体规律的估计。当然，蒙特卡洛方法并不能让不确定性消失。它也无法保证某一次计算结果恰好正确。事实上，蒙特卡洛方法最深刻的地方，恰恰在于它没有假装拥有这种确定性。它允许计算结果存在误差。甚至可以用数学描述误差的规模。它让人类明白：有些问题，我们无法轻易得到精确答案。但这并不意味着我们只能放弃。我们可以通过科学设计的抽样，不断改进自己的认识。有时，最有价值的不是一份声称绝对正确的数字。而是一份知道自己可能错多少的估计。乌拉姆在康复期间思考纸牌概率时，未必能预见这套方法的全部未来。但他所捕捉到的思想，如今已经成为现代计算科学的重要组成部分。上一章，我们发现计算机能够使用确定性算法生成看起来随机的数字。这一章，我们进一步发现：那些数字可以让计算机模拟充满不确定性的世界。从确定性的程序中产生随机样本。再从无数随机样本中发现统计规律。两种看似矛盾的力量，在这里完成了一次结合。人类不能要求每一次偶然都按照自己的意愿发生，但可以让无数次偶然，帮助我们逐渐接近答案。历史档案：事实来源与叙事边界本章主要使用以下历史资料：MacTutor 数学史档案：斯坦尼斯瓦夫·乌拉姆传记。记录其成长经历、1946年患病、纸牌游戏启发与研究生涯。Metropolis 与 Ulam：《The Monte Carlo Method》（1949）。发表于《美国统计协会杂志》第44卷，第335—341页，是蒙特卡洛方法早期重要的公开论文。洛斯阿拉莫斯国家实验室：蒙特卡洛历史文献收藏。收录梅特罗波利斯、Eckhardt 等人关于蒙特卡洛方法起源与发展的回顾。Nicholas Metropolis：Monte Carlo: In the Beginning and Some Great Expectations（1985）。梅特罗波利斯关于方法早期发展、费米的相关探索及计算机应用的回忆。Neutronics Calculation Advances at Los Alamos: Manhattan Project to Monte Carlo（2021）。梳理洛斯阿拉莫斯粒子输运计算与蒙特卡洛发展的历史。本章将乌拉姆的纸牌经历作为有文献支持的当事人回忆；具体康复场景及叙事心理属于文学重构。布丰投针、费米的早期探索、1947年前后的计算研究、1949年的论文和1950年代的后续发展均有史料依据。圆周率估计、简单粒子模型、随机抽样图及 Python 程序为现代教学实验，不是洛斯阿拉莫斯当时的原始计算数据。下一章预告 · 第35章《如何判断一个天文数字是不是素数？》如果给你一个拥有三百位数字的巨大整数，你怎样判断它究竟是不是素数？最直接的方法，是尝试寻找它有没有其他约数。但当数字越来越大时，这条道路可能变得难以承受。二十世纪后期，数学家开始发现，某些数论性质可以帮助我们快速识别合数。更令人惊讶的是，我们甚至可以允许算法使用随机选择，却仍然严格控制它作出错误判断的概率。下一章，我们将从费马小定理出发，经过卡迈克尔数带来的陷阱，走向 Miller–Rabin 素性测试，探索概率算法如何进入现代计算数论与密码学。——《算法的故事》第二季 · 第34章 完上一章下一章