Visualize-ML

💭

"任何考虑用算术手段来产生随机数的人当然都是有原罪的。"

"Anyone who considers arithmetical methods of producing random digits is, of course, in a state of sin."

—— 约翰·冯·诺伊曼 (John von Neumann)
美国籍数学家 | 1903 ~ 1957

🐍 Python工具箱

matplotlib.patches.Circle()

绘制正圆

matplotlib.pyplot.semilogx()

横轴设置为对数坐标

numpy.random.uniform()

均匀分布随机数

numpy.random.normal()

正态分布随机数

numpy.random.binomial()

二项分布随机数

scipy.integrate.quad()

求解定积分值

📋 本章内容概览

基础应用

• 估算平方根
• 估算积分
• 估算体积
• 估算圆周率

抽样方法

• 布丰投针
• 接受-拒绝抽样
• 二项分布随机漫步

高级应用

• 高斯分布相加
• 产生特定相关性
• 数值模拟

✦

15.1

蒙特卡洛模拟：基于伪随机数发生器

Monte Carlo Simulation: Based on Pseudo-Random Number Generators

蒙特卡洛模拟 (Monte Carlo simulation)，也称统计模拟方法，
是以概率统计理论为核心的数值计算方法。

蒙特卡洛模拟将提供多种可能的结果以及通过大量随机数据样本得出的每种结果的概率。
冯·诺伊曼 (John von Neumann) 等三名科学家在20世纪40年代发明了蒙特卡洛模拟。
他们以摩纳哥著名的赌城——蒙特卡洛 (Monte Carlo) ——为其命名。

🎲 伪随机数生成器演示

Pseudo-Random Number Generator Demo

分布类型

样本大小

n = 1000

统计量

样本数

1000

均值

NaN

标准差

NaN

最小值

Infinity

最大值

-Infinity

当前分布信息

均匀分布 U(0, 1)：
• numpy.random.uniform()
• numpy.random.random()
• 理论均值：0.5
• 理论标准差：≈0.289

所有值在[0, 1)区间内等概率出现

前20个样本值

NumPy常用随机函数

random.uniform() - 均匀分布

random.normal() - 正态分布

random.exponential() - 指数分布

random.binomial() - 二项分布

random.poisson() - 泊松分布

random.beta() - 贝塔分布

关键概念：

伪随机数 — 由算术方法产生，确定性但看似随机
种子(seed) — 设置种子可以重现相同的随机序列
大数定律 — 样本量越大，样本统计量越接近理论值
蒙特卡洛 — 基于大量随机样本进行数值计算和模拟

📊 表1. NumPy中和随机数有关的常见函数

函数名称	函数介绍
`numpy.random.beta()`	生成指定形状参数的贝塔分布随机数
`numpy.random.binomial()`	返回给定形状的随机二项分布数组
`numpy.random.chisquare()`	生成指定自由度的卡方分布随机数
`numpy.random.choice()`	随机从给定的数组中选择元素
`numpy.random.dirichlet()`	生成指定参数的狄利克雷分布随机数
`numpy.random.exponential()`	生成指定尺度的指数分布随机数
`numpy.random.gamma()`	生成指定形状和尺度的伽马分布随机数
`numpy.random.lognormal()`	生成指定均值和标准差的对数正态分布随机数
`numpy.random.multivariate_normal()`	生成多元正态分布随机数
`numpy.random.normal()`	生成指定均值和标准差的正态分布随机数
`numpy.random.poisson()`	生成指定均值的泊松分布随机数
`numpy.random.power()`	返回给定形状的随机幂律分布数组
`numpy.random.rand()`	返回一个给定形状的随机浮点数数组，值在0到1之间
`numpy.random.randint()`	返回一个给定形状的随机整数数组，值在给定范围之间
`numpy.random.randn()`	返回一个给定形状的随机浮点数数组，值遵循标准正态分布
`numpy.random.random()`	生成[0, 1)之间的随机数
`numpy.random.seed()`	设置随机数生成器的种子，确保随机数生成的可重复性
`numpy.random.shuffle()`	随机打乱给定的数组
`numpy.random.uniform()`	生成指定范围内的均匀分布随机数

💡 本节小结

蒙特卡洛模拟是一种基于随机抽样的数值计算方法，广泛应用于：

物理模拟（粒子运动、热传导等）
金融工程（期权定价、风险评估等）
统计推断（贝叶斯估计、置信区间等）
优化问题（全局优化、模拟退火等）

本章将介绍几个最基本的蒙特卡洛模拟试验，帮助大家理解这一强大工具的原理和应用。

✦

15.2

估算平方根

Estimating Square Root using Monte Carlo

本节用蒙特卡洛模拟估算√2。

🎲 蒙特卡洛估算√2

Monte Carlo Estimation of Square Root of 2

样本大小

n = 1,000

真实值：

√2 = 1.414214

方法1：直接估算（理论方法，需要已知√2）

原理：在[0, 2]范围内生成随机数，统计小于√2的比例

概率 P(x < √2) = √2 / 2
√2 = 2 × P(x < √2)

样本总数

1,000

x < √2

比例

0.0000

估算结果：

√2 ≈ 0.000000

误差：1.414214 (100.000%)

问题：这个方法需要事先知道√2来判断x是否小于√2，形成循环论证！

⭐ 方法2：平方技巧（实用方法，图1所示）

核心思想：由于√2未知，采用平方技巧避免循环论证

在[0, 2]生成随机数x

计算x² = x × x

统计x² < 2的比例

P(x² < 2) = P(x < √2) = √2 / 2

√2 = √[2 × P(x² < 2)]

样本总数

1,000

x² < 2

比例

0.0000

估算结果：

√2 ≈ 0.000000

误差：1.414214 (100.000%)

✅ 优势：不需要事先知道√2！只需判断x²是否小于2（已知）。

📊 图1解释：平方技巧可视化

原始空间：x ∈ [0, 2]

区间长度：2
目标区间：[0, √2]
目标长度：√2
概率：√2 / 2

问题：√2未知，无法判断

平方空间：x² ∈ [0, 4]

区间长度：4（理论）
实际最大：2² = 4
目标区间：[0, 2]
目标长度：2（已知！）

✅ 可以判断x² < 2

📈 收敛性分析

样本量与精度：
随着样本量增加，估算值会逐渐收敛到真实值√2 = 1.414214...

当前样本量：
1,000

估算精度：
100.000% 误差

💡 尝试增加样本量（如50000），观察误差如何减小

🔢 前10个样本示例

序号	x	x²	x < √2?	x² < 2?

核心要点：

平方技巧 — 通过x²变换，将未知值√2转化为已知值2
概率等价 — P(x < √2) = P(x² < 2)，两个事件概率相同
蒙特卡洛核心 — 用频率估计概率，样本量越大越准确
实用性 — 方法2不需要预先知道答案，可以真正用于估算
代码文件 — Bk5_Ch15_01.py 实现了这个估算过程

💡 理论推导

基本思路

如图1所示，为了估算√2，可以在0 ~ 2的范围内产生大量服从均匀分布的随机数。

在0 ~ 2范围内，随机数在0 ~ √2出现的概率为√2 / 2。

直接方法（理论）

√2则可以根据下式估计得到：

√2 = 2 × [n(x < √2) / n_total]

其中 n() 计算频数。

问题：这个方法需要预先知道√2的值才能判断x是否小于√2，形成循环论证！

⭐ 平方技巧（实用方法）

由于√2未知，所以采用图1所示平方技巧，即√2可以根据下式得到：

关键变换：

x ∈ [0, 2] → x² ∈ [0, 4]

P(x < √2) = P(x² < 2)

√2 = √[2 × n(x² < 2) / n_total]

✅ 优势：只需要判断x²是否小于2（已知值），不需要预先知道√2！

📊 图1. 估算√2的几何解释

原始空间：x = √2

x ∈ [0, 2]

0 ——— 1 ——— √2 ——— 2

目标：找到√2的位置

问题：√2未知，无法判断x < √2

平方空间：x² = 2

x² ∈ [0, 4]（理论）
实际：x² ∈ [0, 2²] = [0, 4]

0 ——— 1 ——— 2 ——— 3 ——— 4

已知：2是固定的阈值

✅ 可以判断x² < 2（2是已知的）

数学等价性

x < √2 ⟺ x² < (√2)² ⟺ x² < 2

两个条件完全等价，但右边的条件可以实际判断！

💻 Python代码

代码文件 Bk5_Ch15_01.py 估算 √2。

🎯 15.2 关键要点

蒙特卡洛方法

• 用频率估计概率
• 基于大量随机样本
• 样本量↑ → 精度↑
• 简单但强大

平方技巧

• 避免循环论证
• x² 变换巧妙
• 已知值判断未知值
• 实用且可行

💡 核心洞察

平方技巧的本质：通过数学变换（x → x²），
将一个需要已知答案才能判断的问题（x < √2？），
转化为一个可以直接判断的问题（x² < 2？）。
这就是蒙特卡洛方法的智慧所在！

✦

15.3

估算积分

Estimating Integrals using Monte Carlo

本节给出的例子用蒙特卡洛模拟方法估算积分。

📊 蒙特卡洛估算积分

Monte Carlo Integration Estimation

样本大小

n = 1,000

积分区间：

[2, 10]

函数最大值：

y_max = 12

函数定义（公式3）

f(x) = (x/2) × sin(x) + 8

在区间 [2, 10] 上的定积分

定积分（公式4）

∫₂¹⁰ [(x/2) × sin(x) + 8] dx

计算函数在 [2, 10] 区间内曲线下的面积

蒙特卡洛方法（图2）

原理：

1. 在长方形区域 [2, 10] × [0, 12] 内随机生成点
2. 判断每个点是否在曲线 f(x) 下方
3. 统计曲线下方的点数比例
4. 用比例乘以长方形面积得到积分估计

长方形面积

8 × 12

样本总数

1,000

曲线下方

蓝色点

曲线上方

1,000

红色点

比例：

0.0000

估算公式（公式5）

A ≈ 矩形面积 × (曲线下方点数 / 总点数)

A ≈ 96 × (0 / 1000)

蒙特卡洛估算结果：

A ≈ 0.0000

精度分析

蒙特卡洛估算

0.0000

数值积分（真值）

67.0526

误差

67.0526

(100.00%)

💡 增加样本数可以提高估算精度，观察误差如何减小

📊 图2说明

如图2所示，在 [2, 10] 区间中，函数 f(x) 的最大值为 12。

在横轴取值从 2 ~ 10，纵轴取值从 0 ~ 12 的长方形空间里，
产生满足均匀分布的 1000 个数据点。

蓝色点 (-)： 在曲线之下

红色点 (×)： 在曲线之上

图2中整个长方形的面积为 96，
定积分对应曲线之下的面积 A，可以通过公式(5)估算得到。

🔢 前10个样本示例

序号	x	y	f(x)	位置

核心要点：

几何意义 — 定积分 = 曲线下方的面积
蒙特卡洛思想 — 用随机点的比例估计面积比例
长方形框定 — 需要知道函数的取值范围（0 ~ 12）
点的分类 — 判断 y ≤ f(x) 确定点是否在曲线下方
精度控制 — 样本量越大，估算越准确
代码文件 — Bk5_Ch15_02.py 实现了这个估算过程

📐 函数定义

给出如下函数 f(x)（公式3）：

f(x) = (x/2) × sin(x) + 8

这是一个周期性振荡的函数，基线在 y = 8 附近

计算 f(x) 在 [2, 10] 区间内定积分（公式4）：

∫₂¹⁰ [(x/2) × sin(x) + 8] dx

计算从 x = 2 到 x = 10 的曲线下方面积

🎲 蒙特卡洛方法原理

步骤1：确定长方形区域

如图2所示，在 [2, 10] 区间中，函数 f(x) 的最大值为 12。

横轴范围： x ∈ [2, 10]
纵轴范围： y ∈ [0, 12]

长方形宽度： 10 - 2 = 8
长方形高度： 12 - 0 = 12
长方形面积： 8 × 12 = 96

步骤2：随机生成点

在横轴取值从 2 ~ 10，纵轴取值从 0 ~ 12 的长方形空间里，
产生满足均匀分布的大量数据点（如1000个）。

步骤3：判断点的位置

对于每个随机点 (x, y)，计算 f(x) 并判断：

蓝色点 (-)
如果 y ≤ f(x)，点在曲线之下

红色点 (×)
如果 y > f(x)，点在曲线之上

步骤4：计算面积（公式5）

图2中整个长方形的面积为 96，
定积分对应曲线之下的面积 A，可以通过下式估算得到：

A ≈ 96 × [n(below f(x)) / 1000]

= 长方形面积 × (曲线下方点数 / 总点数)

n(below f(x)) 为 1000 个数据点中位于 f(x) 曲线之下的数量。

📊 图2. 用蒙特卡洛模拟法求积分

图2展示了蒙特卡洛方法估算积分的可视化过程：

坐标系设置

• 横轴：x ∈ [2, 10]
• 纵轴：y ∈ [0, 12]
• 显示函数曲线 f(x)

随机点分布

• 蓝色点：曲线下方
• 红色点：曲线上方
• 均匀分布在矩形内

关键观察

蓝色点的比例 ≈ 曲线下方面积 / 长方形面积
因此：积分值 ≈ 长方形面积 × 蓝色点比例

💡

数学洞察

蒙特卡洛积分的本质是用面积比估计积分：

定积分 / 长方形面积 ≈ 曲线下方点数 / 总点数

这种方法特别适用于：

难以用解析方法求解的积分
高维积分（多重积分）
复杂区域上的积分

💻 Python代码

代码文件 Bk5_Ch15_02.py 估算积分。

🎯 15.3 关键要点

蒙特卡洛积分法

• 几何意义：定积分 = 面积
• 随机采样：均匀分布的点
• 比例估计：点数比 ≈ 面积比
• 简单有效：易于实现

方法优势

• 不需要求原函数
• 适用于复杂函数
• 可扩展到高维
• 误差可控

💡 核心思想

用随机点的空间分布比例来估计面积比例，
进而计算定积分。这是概率思想在数值计算中的巧妙应用！
样本量越大，估算越准确（大数定律）。

✦

15.4

估算体积

Estimating Volume using Monte Carlo

本节用蒙特卡洛模拟估算空间体积大小。

🎲 蒙特卡洛估算3D体积

Monte Carlo 3D Volume Estimation

样本大小

n = 5,000

空间范围：

x ∈ [-1, 1]
y ∈ [-1, 1]
z ∈ [0, 2]

立方体体积：

二次曲面解析式（公式6）

z = 2 - x² - y²

这是一个开口向下的抛物面

双重定积分（公式7）

∫₋₁¹ ∫₋₁¹ (2 - x² - y²) dx dy

计算曲面与 z = 0 平面之间的体积

蒙特卡洛方法（图3）

步骤：

1. 在立方体空间 [-1,1] × [-1,1] × [0,2] 内随机生成点
2. 对每个点 (x, y, z)，计算曲面高度 z_surface = 2 - x² - y²
3. 判断 z ≤ z_surface（点在曲面下方）
4. 统计曲面下方点的比例
5. 用比例乘以立方体体积得到体积估计

立方体体积

样本总数

5,000

比例

0.0000

曲面下方

蓝色点

曲面上方

5,000

红色点

体积估算结果

V ≈ 立方体体积 × (曲面下方点数 / 总点数)

V ≈ 8 × (0 / 5000)

蒙特卡洛估算体积：

V ≈ 0.0000

精度分析

蒙特卡洛估算

0.0000

数值积分（真值）

5.3336

误差

5.3336

(100.00%)

💡 图3(b)中估算体积约为 5.3712，接近真实值

📊 图3说明

图3(a)：二次曲面

• 显示曲面 z = 2 - x² - y²
• x, y 范围：[-1, 1]
• 蓝色平面：z = 0
• 目标：计算中间体积

图3(b)：随机点分布

• 蓝色点：曲面下方
• 红色点：曲面上方
• 均匀分布在立方体内
• 估算体积 ≈ 5.3712

空间设置

当 x 和 y 均在 [-1, 1] 范围内时，
用蒙特卡洛模拟估算图3(a)曲面和 z = 0 平面（蓝色）构造的空间体积。
整个立方体空间体积为 8，在这个空间均匀产生大量随机点（如5000个）。

方法对比

类似上一节（15.3估算积分），
根据随机点的比例，可以估算公式(7)的双重定积分。
区别：从2D积分扩展到3D体积估算！

🔢 前10个样本示例

序号	x	y	z	z_surface	位置

核心要点：

3D扩展 — 从2D积分扩展到3D体积估算
双重积分 — ∬(2 - x² - y²) dx dy 表示体积
立方体框定 — 需要确定包含曲面的立方体空间
点的判断 — z ≤ 2 - x² - y² 确定点在曲面下方
比例估算 — 体积 ≈ 立方体体积 × 下方点比例
代码文件 — Bk5_Ch15_03.py 估算本节体积

🎨 二次曲面定义

图3(a)所示二次曲面解析式（公式6）：

z = 2 - x² - y²

这是一个开口向下的抛物面（paraboloid）

顶点在 (0, 0, 2)，沿z轴负方向延伸

空间范围

当 x 和 y 均在 [-1, 1] 范围内时，
估算图3(a)曲面和 z = 0 平面（蓝色）构造的空间体积。

x 范围： [-1, 1]
y 范围： [-1, 1]
z 范围： [0, 2]

立方体体积：
2 × 2 × 2 = 8

∬ 双重定积分（公式7）

体积的数学表达

这个体积相当于如下双重定积分：

V = ∫₋₁¹ ∫₋₁¹ (2 - x² - y²) dx dy

对 x 和 y 在 [-1, 1] 区间内进行二重积分

积分的几何意义

理解方式1：

对于每个 (x, y) 点，函数值 2 - x² - y² 表示该点处的高度。
对所有点的高度求和（积分），就得到总体积。

理解方式2：

想象把区域 [-1,1] × [-1,1] 分成无数小方块，
每个小方块上方的"柱子"高度是 2 - x² - y²，
所有柱子的体积之和就是总体积。

🎲 蒙特卡洛3D方法

方法步骤

确定包围盒：
立方体空间 [-1, 1] × [-1, 1] × [0, 2]，体积 = 8
生成随机点：
在立方体内均匀产生大量随机点（如5000个）
判断位置：
对每个点 (x, y, z)，计算 z_surface = 2 - x² - y²
如果 z ≤ z_surface，点在曲面下方
统计比例：
计算曲面下方点数占总点数的比例
估算体积：
体积 ≈ 立方体体积 × 曲面下方点比例

点的分类（图3b）

蓝色点

z ≤ 2 - x² - y²
点在曲面之下
这些点对应我们要计算的体积

红色点

z > 2 - x² - y²
点在曲面之上
这些点在立方体内但不在目标体积内

📊 图3. 利用蒙特卡洛估算体积

图3(a)：二次曲面

曲面展示：
• 3D图像显示抛物面 z = 2 - x² - y²
• 曲面颜色可能渐变显示高度
• 底部蓝色平面表示 z = 0

几何特征：
• 顶点：(0, 0, 2)
• 对称性：关于 z 轴对称
• 开口：向下（负 z 方向）

图3(b)：随机点分布

可视化：
• 3D散点图显示随机点
• 蓝色点密集在曲面下方
• 红色点稀疏在曲面上方

估算结果：
Estimated volume = 5.3712
（图中标注的估算值）

从2D到3D的扩展

15.2 估算√2：
1D → 线段
判断 x < √2

15.3 估算积分：
2D → 面积
判断 y < f(x)

15.4 估算体积：
3D → 体积
判断 z < f(x,y)

统一思想：用随机点的比例估计几何量的比例！

💻 Python代码

Bk5_Ch15_03.py 估算本节体积。

🎯 15.4 关键要点

3D蒙特卡洛

• 双重积分 → 体积
• 立方体包围盒
• 3D随机点生成
• 曲面位置判断

方法扩展

• 从2D到3D自然扩展
• 原理完全一致
• 可推广到更高维
• 计算复杂度可控

💡 核心洞察

蒙特卡洛方法的强大之处在于维度扩展性：
无论是1D、2D还是3D，甚至更高维，
核心思想都是一样的——用随机采样的比例来估计几何量。
这使得蒙特卡洛成为处理高维积分问题的有力工具！

✦

15.7

接受-拒绝抽样法

Accept-Reject Sampling Method

随机变量X的概率密度函数为 f_X(x)，
但是 f_X(x) 不可以直接抽样。
我们可以采用本节介绍的接受-拒绝抽样法 (accept-reject sampling method)。

✂️ 接受-拒绝抽样法

Accept-Reject Sampling Method

样本数量

n = 2,000

辅助分布：

Uniform[-4, 4] × [0, 0.5]

接受率：

NaN%

方法概述

随机变量X的概率密度函数为 f_X(x)，但是 f_X(x) 不可以直接抽样。
我们可以采用接受-拒绝抽样法。

适用场景：
概率密度函数复杂、不能直接抽样的情况

基本思想

1. 生成辅助分布

生成一个辅助分布 (proposal distribution)，
并利用这个分布来生成随机数。

2. 计算接受率

计算目标概率分布在该点处的概率密度，
并将其除以辅助分布在该点处的概率密度，
得到接受率 (acceptance ratio)。

3. 判断接受或拒绝

随机生成一个介于0和1之间的均匀分布的随机数，
如果这个随机数小于接受率，则接受这个样本，
否则拒绝。重复此过程，直到生成足够多的样本。

图10-13：可视化过程

图10. 目标概率密度函数

给定随机变量X的概率密度函数 f_X(x)。
显然没有"现成"的随机数发生器能够直接生成满足 f_X(x) 的随机数。

图11. 生成连续均匀分布随机数

首先生成连续均匀分布随机数。
在矩形区域 [x_min, x_max] × [0, y_max] 内均匀采样。

图12. "剪裁"连续均匀分布随机数

简单来说，接受-拒绝抽样法就是在图11中"剪裁"得到形似图10的部分，
并"接受"这些随机数。

图13. "接受" vs "拒绝"部分随机数

用直方图可视化"拒绝"和"接受"部分的随机数。
浅蓝色部分矩形构成形状形似图10中的 f_X(x)。

抽样结果统计

总样本数

2,000

接受样本

浅蓝色

拒绝样本

深色

接受率：

NaN%

= 0 / 2000 × 100%

优点与缺点

✓ 优点

• 简单易用
算法原理直观，容易理解和实现

• 适用范围广
可以应用于各种不同的概率分布

✗ 缺点

• 效率可能较低
需要进行接受和拒绝的判断

• 辅助分布选择重要
对样本生成的效率和精度影响大

📊 可视化说明

接受的样本（绿色）：
点在目标PDF曲线下方，y ≤ f_X(x)

拒绝的样本（红色）：
点在目标PDF曲线上方，y > f_X(x)

观察：
接受的样本在x轴上的分布会逐渐逼近目标概率密度函数 f_X(x)。
图13的直方图浅蓝色部分形状会形似图10的目标PDF。

🔢 前10个样本示例

序号	x	y	f_X(x)	判断

📖 后续应用

我们将在本书第22章用到接受-拒绝抽样法。

核心要点：

问题场景 — 目标PDF复杂，无法直接抽样
辅助分布 — 生成易于采样的均匀分布
剪裁思想 — 从均匀分布中"剪裁"出目标分布
接受条件 — y ≤ f_X(x)则接受，否则拒绝
效率考量 — 接受率越高，效率越高
图形直观 — 图10-13完整展示了整个过程

❓

问题场景

随机变量X的概率密度函数为 f_X(x)，
但是 f_X(x) 不可以直接抽样。

也就是说，不能直接产生满足 f_X(x) 的随机数。

解决方案：
接受-拒绝抽样法适合于概率密度函数复杂、不能直接抽样的情况。

💡 基本思想

步骤1：生成辅助分布

接受-拒绝抽样法的基本思想是，
生成一个辅助分布 (proposal distribution)，
并利用这个分布来生成随机数。

步骤2：计算接受率

然后，计算目标概率分布在该点处的概率密度，
并将其除以辅助分布在该点处的概率密度，
得到接受率 (acceptance ratio)。

步骤3：判断接受或拒绝

随机生成一个介于0和1之间的均匀分布的随机数，
如果这个随机数小于接受率，则接受这个样本，
否则拒绝。
重复此过程，直到生成足够多的样本。

📊 图10-13：可视化过程

图10. 随机变量X的概率密度函数

给定如图10所示的随机变量X的概率密度函数为 f_X(x)。
显然没有"现成"的随机数发生器能够直接生成满足 f_X(x) 的随机数。

图11. 生成连续均匀分布随机数

我们首先生成如图11所示的连续均匀分布随机数。
在矩形区域内均匀采样：x ∈ [x_min, x_max]，y ∈ [0, y_max]

图12. "剪裁"连续均匀分布随机数

简单来说，如图12所示，接受-拒绝抽样法就是在图11中
"剪裁"得到形似图10的部分，并"接受"这些随机数。

剪裁规则：
对于每个随机点 (x, y)：
• 如果 y ≤ f_X(x)，则接受该样本
• 如果 y > f_X(x)，则拒绝该样本

图13. "接受" vs "拒绝"部分随机数

图13所示为用直方图可视化"拒绝"和"接受"部分的随机数。
大家很容易发现，图中浅蓝色部分矩形构成形状形似图10中的 f_X(x)。

这正是接受-拒绝抽样法的精髓：通过"剪裁"获得目标分布！

⚖️ 优点与缺点

✓优点

简单易用

算法原理直观，容易理解和实现

适用范围广

可以应用于各种不同的概率分布

✗缺点

样本生成效率可能较低

因为需要进行接受和拒绝的判断

辅助分布选择重要

在实践中，辅助分布的选择对于样本生成的效率和精度非常重要

🎨

几何直观理解

接受-拒绝抽样法的几何意义非常直观：

1. 均匀撒点：在矩形区域内均匀随机撒点
2. 曲线判断：判断每个点是否在目标PDF曲线下方
3. 接受下方：保留曲线下方的点，丢弃上方的点
4. 形成分布：接受的点在x轴上的投影形成目标分布

这就像用"筛子"从均匀分布中筛选出目标分布的样本！

📖 后续应用

我们将在本书第22章用到接受-拒绝抽样法。

🎯 15.7 关键要点

方法核心

• 问题：目标PDF无法直接抽样
• 思想：从均匀分布"剪裁"
• 判断：y ≤ f_X(x)则接受
• 结果：形成目标分布

实践要点

• 辅助分布选择很关键
• 接受率影响效率
• 适用范围广泛
• 原理简单直观

💡 核心洞察

接受-拒绝抽样法是蒙特卡洛方法的重要技术：
• 将复杂问题转化为简单的几何判断
• 用"剪裁"思想从简单分布获得复杂分布
• 图10-13完美展示了从概念到实现的全过程
• 在后续章节（第22章）会有更深入的应用！

✦

15.9

两个服从高斯分布的随机变量相加

Sum of Two Gaussian Random Variables

X₁ 和 X₂ 分别服从正态分布，
它们的和也服从正态分布。

📊 两个高斯分布随机变量相加

Sum of Two Gaussian Random Variables

X₁参数

μ₁ = -1

σ₁ = 1

X₂参数

μ₂ = 1

σ₂ = 1

🎯 15.9 关键要点

正态分布性质

• 两正态分布相加仍为正态分布
• 均值：μY = μ₁ + μ₂
• 方差：σY² = σ₁² + σ₂² + 2ρσ₁σ₂
• 标准差：σY = √(σY²)

产生满足特定相关性的随机数

Generate Correlated Random Numbers via Cholesky Decomposition

利用Cholesky分解，
我们可以产生满足特定相关性的随机数。

🔢 产生满足特定相关性的随机数

Cholesky Decomposition for Correlated Random Numbers

期望值 μ

μ₁ = 2

μ₂ = 4

标准差 σ

σ₁ = 2.00

σ₂ = 1.41

与《矩阵力量》的联系（图23）

如图23所示，我们在《矩阵力量》第14章中学过如何完成
"单位圆（缩放）→ 正椭圆（剪切）→ 旋转椭圆"几何变换。

经过本书上一板块的学习，大家已经清楚：
• 单位圆代表 N(0, I)
• 旋转椭圆代表 N(0, Σ)
• 再经过平移，我们就可以到 N(μ, Σ)

图23用到的数学工具是LDL分解。
实际上，利用Cholesky分解我们可以通过一次矩阵乘法
便完成"缩放 + 剪切"。
这便是利用Cholesky分解结果产生满足特定相关性随机数的技术路线。

📊 图24. 矩阵运算示意图

步骤1：生成标准正态随机数矩阵

如图24所示，首先生成满足 N(0, I_D×D) 的随机数矩阵 Z。

Z_n×D ~ N(0, I_D×D)

步骤2：Cholesky分解协方差矩阵

然后，对协方差矩阵 Σ 进行Cholesky分解，
得到下三角矩阵 R^T 和上三角矩阵 R 乘积：

Σ = R^TR （公式21）

注意：要求Σ为正定；否则，不能进行Cholesky分解。

步骤3：矩阵乘法与平移

矩阵R中含有图23中"剪切"、"缩放"两个成分。

X = ZR + E(X) （公式22）

Z服从N(0, I_D×D)，
经过如下运算得到的多元随机数X服从N(E(X), Σ_D×D)

🧮 公式23-25：协方差矩阵验证

向量定义（公式23）

X = [x₁ x₂ ... x_D]
Z = [z₁ z₂ ... z_D]
E(X) = [μ₁ μ₂ ... μ_D]

Z的协方差矩阵（公式24-25）

Σ_Z = (1/(n-1)) Z^TZ = I

因为Z的每一列都是独立的标准正态分布，
所以协方差矩阵是单位矩阵I。

X的协方差矩阵推导

Σ_X = (1/(n-1)) (X - E(X))^T(X - E(X))
= (1/(n-1)) (ZR)^T(ZR)
= (1/(n-1)) R^TZ^TZR
= R^T · I · R
= R^TR = Σ ✓

验证成功！X的协方差矩阵确实等于Σ

📝 二维情况（公式26-29）

目标分布（公式26）

二维随机变量χ满足如下二维高斯分布：

χ = [X₁, X₂]^T ~ N(μ, Σ)

μ = [μ₁, μ₂]^T

Σ = [σ₁², ρσ₁σ₂]
[ρσ₁σ₂, σ₂²]

标准正态分布（公式27）

而Z₁和Z₂服从标准正态分布，且不相关：

ς = [Z₁, Z₂]^T ~ N(0, I₂ₓ₂)

Σ_Z = [1, 0]
[0, 1]

Cholesky分解（公式28）

对协方差矩阵Σ进行Cholesky分解：

L = [σ₁, 0]
[ρσ₂, σ₂√(1-ρ²)]

变换关系（公式29）

χ与ς的关系为：

χ = Lς + μ

思考：为什么公式22采用上三角矩阵R，
而公式29采用下三角矩阵L？
（答：取决于矩阵乘法的顺序，ZR或Lς）

展开形式（公式32-33）

D=2时，展开公式29得到：

X₁ = σ₁·Z₁ + μ₁
X₂ = ρσ₂·Z₁ + σ₂√(1-ρ²)·Z₂ + μ₂

这就是从标准正态分布生成相关正态分布的具体公式！

📊

图25. 具体示例

下面给出一个具体示例：

图25 (a)：

[Z₁, Z₂]^T ~ N([0,0], I)

二维随机数满足标准正态分布

图25 (b)：

[X₁, X₂]^T ~ N([2,4], Σ)

Σ = [4, 2]
[2, 2]

满足特定质心和协方差要求

将服从IID二维标准正态分布随机数
转化为满足特定质心和协方差要求的随机数

🔢 多维随机数（图26）

图26所示为采用多元高斯分布随机数发生器生成的随机数。
这组随机数的均值、协方差矩阵和鸢尾花数据相同。
请大家利用本节前文介绍的技术原理：
• 首先生成满足 N(0, I_D×D) 的随机数矩阵Z
• 然后再生成满足 N(E(X), Σ_D×D) 的随机数

💻 Python代码

Bk5_Ch15_08.py 生成图25。代码中用到了Cholesky分解。

🎯 15.10 关键要点

Cholesky分解

• Σ = R^TR 或 Σ = LL^T
• 要求Σ正定
• 一次矩阵乘法完成缩放+剪切
• 下三角L或上三角R

变换公式

• X = ZR + E(X)
• Z ~ N(0, I) → X ~ N(μ, Σ)
• 图24完整流程
• 图25二维示例

💡 核心洞察

Cholesky分解是产生相关随机数的核心技术：
• 连接了《矩阵力量》的几何变换与蒙特卡洛模拟
• 从简单的标准正态分布生成复杂的相关分布
• 协方差矩阵的分解体现了"缩放+剪切"的几何意义
• 图26展示了多维扩展的可能性（如鸢尾花数据）
• 这是多元统计分析和机器学习中的基础技术！

🐍 Python工具箱

📋 本章内容概览

蒙特卡洛模拟：基于伪随机数发生器

🎲 伪随机数生成器演示

分布类型

样本大小

📊 表1. NumPy中和随机数有关的常见函数

💡 本节小结

估算平方根

🎲 蒙特卡洛估算√2

样本大小

💡 理论推导

📊 图1. 估算√2的几何解释

🎯 15.2 关键要点

估算积分

📊 蒙特卡洛估算积分

样本大小

📐 函数定义

🎲 蒙特卡洛方法原理

📊 图2. 用蒙特卡洛模拟法求积分

数学洞察

🎯 15.3 关键要点

估算体积

🎲 蒙特卡洛估算3D体积

样本大小

🎨 二次曲面定义

∬ 双重定积分（公式7）

🎲 蒙特卡洛3D方法

📊 图3. 利用蒙特卡洛估算体积

🎯 15.4 关键要点

接受-拒绝抽样法

✂️ 接受-拒绝抽样法

样本数量

问题场景

💡 基本思想

📊 图10-13：可视化过程

⚖️ 优点与缺点

几何直观理解

🎯 15.7 关键要点

两个服从高斯分布的随机变量相加

📊 两个高斯分布随机变量相加

X₁参数

X₂参数

相关系数

🎯 15.9 关键要点

产生满足特定相关性的随机数

🔢 产生满足特定相关性的随机数

期望值 μ

标准差 σ

相关系数 ρ

与《矩阵力量》的联系（图23）

📊 图24. 矩阵运算示意图

🧮 公式23-25：协方差矩阵验证

📝 二维情况（公式26-29）

图25. 具体示例

🎯 15.10 关键要点