【统计方法】蒙特卡洛

Monte Carlo方法详解

什么是Monte Carlo方法？

Monte Carlo方法是一类依赖随机采样的计算技术，广泛应用于各种难以解析求解的问题。

• 主要特点是依赖随机采样来近似解决问题，提供的是近似解而非精确解。
• 常用于解析解不存在或难以实现的场景。
• 有时也称为随机模拟（stochastic simulation）。

Monte Carlo方法的起源

• 由波兰数学家斯坦尼斯瓦夫·乌拉姆（Stanislaw Ulam）在20世纪40年代末发明。
• 灵感来源于他病中思考纸牌游戏“Canfield Solitaire”成功概率的估计，
• 通过重复实验而非纯组合计算来估计概率。
• “Monte Carlo”是项目代号，来源于摩纳哥著名赌场，寓意随机性和赌博。

Monte Carlo估计$\pi$

利用随机点在单位正方形中落入单位圆的比例来估计$\pi$。

• 单位圆面积为$\pi r^2=\pi$（$r=1$），单位正方形面积为$4r^2=4$。
• 随机生成$N$个点，计算落入单位圆内的点数$R$。
• 估计$\pi \approx 4\frac{R}{N}$。

统计学原理

将目标参数写成期望形式：
$$\pi =\mathbb{E}[X]={\int }_{A}tf(t)dt\approx \frac{1}{N}\sum _{i=1}^N{X}_i$$
其中$X_i$是采样的随机变量，$f(\cdot )$为概率密度函数，$A$为定义域。
根据大数定律，随着样本数量增加，估计值趋近真实值。

• 这表示在区域$A$上，所有可能值$t$按概率$f(t)$加权求和（积分）。
• 期望值是随机变量的“加权平均值”

这是一个通用的期望估计公式，表示从分布$f$采样$X_i$后，$X_i$的平均值近似期望。

• 但是，估计$\pi$的Monte Carlo方法中，$X_i$本身并不是直接的点坐标，而是一个指示变量：

$$X_i=\{\begin{array}{ll}1,& \text{点落在圆内}\\ 0,& \text{点落在圆外}\end{array}$$

• 于是，这时期望$\mathbb{E}[X]$表示“点落在圆内的概率”，其实是：

$$\mathbb{E}[X]=P(\text{点落在圆内})=\frac{\pi }{4}$$

因为采样是在正方形内，面积是4，而圆面积是$\pi$。

一般期望的近似

若有函数$g(X)$，同样可估计其期望：
$$\mathbb{E}[g(X)]={\int }_{A}g(t)f(t)dt\approx \frac{1}{N}\sum _{i=1}^{N}g(X_i)$$
假设$X_i$从分布$f$中独立抽样。

• 假设我们有一组样本$X_i$服从$f$分布。
• 计算$\mathbb{E}[X^2]$，即$g(x)=x^2$的期望。
• 样本均值计算为

$$\frac{1}{n}\sum _{i=1}^n{X}_i^2$$

• 即用样本平方求平均，作为$\mathbb{E}[X^2]$的估计。

关键理解点

1. 期望是积分，但难计算时用随机采样替代积分。
2. 蒙特卡洛方法利用大量独立样本，样本均值近似期望。
3. 函数$g$可以灵活选择，实现对不同函数期望的估计。

Monte Carlo积分示例

目标是估计积分：
$${\int }_0^1{e}^{-x^2/2}dx\approx \frac{1}{N}\sum _{i=1}^{N}g(X_i)$$
其中$g(x)=e^{-x^2/2}$，$X_i$服从$U(0,1)$均匀分布。

R代码实现：

f <- function(x) {
  exp(-x^2/2)
}
integrate(f, lower=0, upper=1)

set.seed(66)
x <- runif(1000, 0, 1)
mean(f(x))

数值结果与精确积分非常接近。

模拟随机变量

• 已知如何模拟均匀分布$U$的随机变量，可用R语言中的runif函数实现。
• 但如何模拟任意概率分布的随机变量？
• 常用两种方法：

1. 逆变换法（Inverse-transform method）
2. 拒绝采样法（Acceptance-rejection method）

逆变换法

• 设随机变量$X$的累计分布函数（CDF）为
  $$F(x)={\int }_{-\infty }^{x}f(t)dt$$

• 记$U\sim \text{Uniform}(0,1)$。

• 若$F^{-1}$存在，则
  $$X=F^{-1}(U)$$
  即通过对均匀变量取CDF的反函数生成目标分布变量。

  逆变换法的核心思想

  • 你知道 $F(X)$ 的取值范围是 $[0,1]$，因为累计概率最大是1，最小是0。
  • 假设 $U$ 是一个在 $[0,1]$ 上均匀分布的随机变量（比如用 runif 函数生成）。
  • 逆变换法告诉我们：

  $X=F^{-1}(U)$

  就是说，如果你用均匀分布 $U$ 取样，再用 $F$ 的反函数把 $U$ 转换回来，就得到了分布为 $X$ 的样本。

  模拟任意分布的随机变量

  在实际问题中，很多随机现象并不是均匀分布，而是其他复杂的分布，比如指数分布、正态分布、卡方分布等等。

  • 逆变换法给出了一种通用的方法：
    • 先在 $[0,1]$ 区间内产生均匀分布随机数（用计算机很容易实现）。
    • 再利用目标分布的累计分布函数（CDF）的反函数，把均匀随机数“转换”为目标分布的随机数。

  这让计算机可以“模拟”各种复杂的随机事件。

指数分布示例

• 指数分布概率密度函数（PDF）：
  $$f(x)=\{\begin{array}{ll}\lambda e^{-\lambda x}& x>0\\ 0& x\le 0\end{array}$$

• 累积分布函数（CDF）：
  $$F(x)=\{\begin{array}{ll}1-e^{-\lambda x}& x>0\\ 0& x\le 0\end{array}$$

• 反函数为：
  $$F^{-1}(p)=-\frac{\log (1-p)}{\lambda },0<p<1$$

  F(x) = U

  $1-e^{-\lambda x}$ = U

  $1-U=e^{-\lambda x}$

  $log(1-U)=-\lambda X$

• 生成指数分布随机变量步骤：

  1. 采样$U\sim \text{Uniform}(0,1)$
  2. 计算$X=-\frac{\log (1-U)}{\lambda }$

拒绝采样法（Acceptance-Rejection method）

• 当逆变换法难以使用时，拒绝采样法更通用。

• 给定两个概率密度$f(x)$和$g(x)$，且满足
  $$f(x)\le Mg(x),\forall x$$

• $g(x)$易于采样。

• 采样过程：

  1. 从$g(x)$采样变量$X$。
  2. 以概率$\frac{f(X)}{Mg(X)}$接受$X$，否则拒绝。
  3. 重复步骤1和2直到获得足够样本。

  1. 概率接受的过程

  • 你采样了一个 X∼g(x)
  • 生成一个均匀随机数 U∼Uniform(0,1)
  • 如果 $\le \frac{f(X)}{Mg(X)}$，接受 X。
  • 否则，拒绝 X，重新采样。

拒绝采样示意图

• $f(x)$是目标分布，通常复杂难采样。
• $Mg(x)$为包络函数，形状包住$f(x)$，常为简单分布的放缩。
• 采样效率依赖于$M$和包络函数对目标分布的贴合程度。

拒绝采样问题

• $g(x)$需与$f(x)$形状相似，否则采样效率低。
• 高维数据中采样效率急剧降低。
• 当维度约为10时，需要采样大量候选样本，效率很低。

总结

• Monte Carlo方法通过随机采样解决积分和概率计算问题。
• 逆变换法适合CDF可逆的分布。
• 拒绝采样法更通用，但效率依赖包络函数。
• Monte Carlo方法在科学计算、机器学习、统计模拟中广泛应用。

Monte Carlo模拟示例：Pop Mart盲盒问题

每个Pop Mart Monster盲盒售价为22美元。
完整系列共有10个独特的模型。
每个盒子随机包含一个模型。

问题：你需要购买多少盒子，预计要花多少钱，才能集齐这10个模型？

解析解（蒙特卡洛模拟前）

期望购买盒子数量计算公式为：
$$E[T]=n\left(\frac{1}{1}+\frac{1}{2}+\cdots +\frac{1}{n}\right)\approx n\log n+0.5772n+0.5$$

• $T$表示需要购买的盒子数量，
• $n$表示总模型数量。

对于$n=10$，计算得到：
$$E[T]\approx 30$$
即平均需要购买约30个盒子。

总预计花费为：
$$30\times 22=660\text{美元}$$

Monte Carlo模拟代码示例（R语言）

n.shops.sim <- c()
for(i in 1:1000) {
  collected <- c()
  n.shops <- 0
  while(length(collected) < 10) {
    newitem <- sample(10, 1)
    collected <- union(collected, newitem)
    n.shops <- n.shops + 1
  }
  n.shops.sim <- c(n.shops.sim, n.shops)
}
mean.sim <- mean(n.shops.sim) # 模拟平均购买盒子数

ggplot(data.frame(x = n.shops.sim)) +
  theme_minimal() +
  geom_histogram(aes(x = x, y = after_stat(density)),
                 colour = "black", fill = "white") +
  labs(x = "Number of Blind Box", y = "Density") +
  geom_vline(xintercept = mean.sim, linetype = "dotted", colour = "red")

初始化一个空列表 n_shops_sim 用于存储每次模拟需要购买的盒子数量

重复进行 1000 次模拟：
初始化一个空集合 collected，用来存放已收集到的独特盒子编号
初始化计数器 n_shops = 0，记录本次模拟购买盒子的数量

当 collected 中盒子数量小于 10 时，循环执行：
随机从 1 到 10 中选择一个盒子编号 new_item
将 new_item 加入到 collected 集合（去重）
购买次数计数器 n_shops 增加 1

将本次模拟的购买次数 n_shops 添加到 n_shops_sim 列表中

计算 n_shops_sim 列表的平均值 mean_sim，作为模拟的期望购买盒子数

结论

• Monte Carlo模拟为复杂概率问题提供数值解法，适合实际随机过程的模拟与评估。
• 模拟结果验证了解析解的准确性和合理性。
• 盲盒集齐问题的模型与实际生活场景紧密相关，具有很强的应用价值。

olya Urn问题

问题背景

• 你有一个装了1个黑球和1个白球的袋子。
• 每次从袋子中抽一个球出来，抽到什么颜色的球，就往袋子里放回“两个”同色球（即原本的球 + 一个新球）。
• 问经过很多次抽取后，黑球和白球的比例会怎样？

模拟思路

1. 初始袋子有 1黑 + 1白。
2. 进行多次抽取：
   • 计算袋子中黑球和白球的数量，求黑球比例。
   • 按比例随机抽取一个球。
   • 抽中球的颜色，加两个该颜色球回袋子。
3. 记录每次抽取后黑球比例。
4. 反复多次模拟，观察比例的变化轨迹。

意义

• 这是一个“自强化”过程，抽到某色球的概率随该色球数量增加而增加，模拟展示这种动态过程。
• 能帮助理解概率模型、随机过程的演变。

欧式期权定价

问题背景

• 欧式看涨期权：到期时，可以按固定价格买入股票。
• 股票当前价格假设100元，执行价105元，90天后看股票价格如何。
• 如果价格高于105元，期权收益 = 股票价格 - 105元，否则收益为0。

如何估价

1. 假设股票价格服从随机过程（几何布朗运动），有增长率 $\mu$ 和波动率 $\sigma$。
2. 利用蒙特卡洛方法模拟股票未来价格的众多可能路径。
3. 对每条路径，计算期权到期收益 $\max (S_T-K,0)$。
4. 取所有模拟路径收益的平均，得到期权的“合理价格”。

模拟步骤伪代码

• for i in 1 to N simulations:
  • 根据随机过程模拟未来股价 $S_T^{(i)}$。
  • 计算收益 $payof{f}^{(i)}=\max (S_T^{(i)}-K,0)$。
• 期权价格 = 所有收益的平均值。

作用

• 解决经典金融定价问题。
• 当解析解困难时，用蒙特卡洛数值估计价格。

股票价格随机模型

模型描述
股票价格未来变化由以下随机微分方程给出：
$$dS=S(\mu dT+\sigma \sqrt{dT}\mathcal{N})$$

• $\mu$：股票预期增长率
• $\sigma$：股票波动率（风险度量）
• $\mathcal{N}$：标准正态随机变量
• $dT$：时间步长（通常为一天或其他单位）

模拟步骤

• 给定当前股价 $S_0$。

• 逐步模拟未来股价：
  $$S_{t+1}=S_t\times (1+\mu dT+\sigma \sqrt{dT}Z)$$
  其中 $Z\sim N(0,1)$。

• 重复生成多条价格路径，统计终值。

意义

• 量化股票价格未来不确定性。
• 为期权等衍生品定价提供基础。

Markov Chain Monte Carlo (MCMC)

Monte Carlo Simulation回顾

• Monte Carlo模拟通过独立从已知分布生成随机样本，来估计期望或概率。
• 示例：模拟股票价格、掷骰子、估算$\pi$等。
• 缺点：
  • 当目标是稀有事件或小区域时效率低。
  • 维度升高时效率急剧下降（维度灾难）。
  • 采样之间无依赖，缺乏“学习”或反馈机制。

Markov Chain Monte Carlo (MCMC)简介

• MCMC是一类用于从任意复杂分布中采样的蒙特卡洛方法。
• 样本不再相互独立，而是形成一个马尔可夫链，后续状态只依赖当前状态。
• 两种常见算法：
  • Metropolis-Hastings算法及其推广
  • Gibbs采样器

马尔可夫链（Markov Chain）

• 马尔可夫链是满足马尔可夫性质的随机过程。

• 马尔可夫性质：未来状态只依赖当前状态，与过去状态无关。

• 考虑随机变量序列 { X 1 , X 2 , … , X n } \{X_1, X_2, \dots, X_n\} {X1​,X2​,…,Xn​}，满足
  $$P(X_n|X_{n-1},X_{n-2},\dots ,X_1)=P(X_n|X_{n-1})$$

马尔可夫链三个关键组成部分

1. 状态空间：所有可能状态的集合。
2. 转移矩阵：从一个状态转移到另一个状态的概率矩阵。
3. 随机序列：状态的随机演变序列 { X 0 , X 1 , X 2 , …   } \{X_0, X_1, X_2, \dots\} {X0​,X1​,X2​,…}。

马尔可夫状态转移示意——天气模型

• 状态空间： S = { Sunny , Cloudy , Rainy } S = \{\text{Sunny}, \text{Cloudy}, \text{Rainy}\} S={Sunny,Cloudy,Rainy}
• 转移概率矩阵$P$：

$$P=\left(\begin{array}{ccc}0.6& 0.1& 0.3\\ 0.4& 0.5& 0.1\\ 0.2& 0.5& 0.3\end{array}\right)$$

• 行表示当前状态，列表示下一状态。
• 例如，今天晴天，明天下雨的概率是0.3。

马尔可夫链长期行为

• 随着时间推移，状态序列会趋于稳定，达到不变分布（stationary distribution）$\pi$，满足：

$$\pi =\pi P,\sum _i{\pi }_i=1$$

• 不变分布代表系统长期的概率分布。

不变分布存在的条件

• 不可约性（Irreducible）：所有状态之间能相互到达。

  • 对一个马尔可夫链或其转移矩阵$P$，当且仅当以下条件满足时，称其为不可约的：

  对于任意状态$i$和$j$，存在一种方式（路径）能从状态$i$达到状态$j$，且也能从状态$j$达到状态$i$

  • 换句话说，所有状态都是互相连通的。

• 非周期性（Aperiodic）：状态返回不受固定周期限制。

  • 考虑以下转移矩阵：
    $$P=\left(\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ 1& 0& 0\end{array}\right)$$

  • 马尔可夫链状态循环经过：
    $$1\to 2\to 3\to 1\to 2\to 3\to \cdots$$

  • 你只能在经过3、6、9……步后返回起始状态。

  • 因此，返回起始状态步数的最大公约数（gcd）为：
    $$\gcd (3,6,9,\dots )=3$$

  • 这说明所有状态都是周期为3的周期状态（periodic with period 3）。

  • 假设我们修改了转移矩阵为：
    $$P=\left(\begin{array}{ccc}0.1& 0.9& 0\\ 0& 0.5& 0.5\\ 0.4& 0& 0.6\end{array}\right)$$

  • 现在状态可以在第1、2、3、4……步返回。

  • 返回步数的最大公约数：
    $$\gcd (1,2,3,4,\dots )=1$$

  • 因此，这些状态是非周期性的（aperiodic）。

不变马尔可夫链的应用

• 和Monte Carlo方法类似，我们的主要兴趣是模拟序列 $X_1,X_2,\cdots$，这些序列服从（或近似服从）某个概率密度函数 $f(x)$，但直接从 $f(x)$ 采样非常困难。

• 这个问题在贝叶斯统计中非常常见：

  • 频率派统计（Frequentist statistics）：存在固定但未知的真实参数 $\theta$。

  • 贝叶斯统计（Bayesian statistics）：所有模型参数都是随机变量。

• 对于贝叶斯统计，模型参数 $\theta$ 有一个先验分布。

• 给定数据 $\mathbf{D}$，我们想要计算后验分布 $p(\theta |\mathbf{D})$：

  $$p(\theta |\mathbf{D})=\frac{p(\mathbf{D}|\theta )p(\theta )}{\int p(\mathbf{D}|{\theta }^{\prime })p({\theta }^{\prime })d{\theta }^{\prime }}$$

• 通常很难精确计算后验分布 $p(\theta |\mathbf{D})$。

• 解决方法是使用MCMC从后验分布 $p(\theta |\mathbf{D})$ 采样。

• 思路是构造一个马尔可夫链，使其平稳分布（stationary distribution）就是我们想要的后验分布 $p(\theta |\mathbf{D})$。

MCMC的核心思想

• 构造一个马尔可夫链，其平稳分布（stationary distribution）正是我们想采样的目标分布。
• 通过模拟该马尔可夫链，生成的样本近似服从目标分布。
• 代表算法有Metropolis-Hastings和Gibbs采样。

Metropolis-Hastings 算法

直观理解

想象你是个政治家，想按选民比例拜访所有城镇。

• 目标：按各城镇选民比例在每个城镇停留的时间成比例。

工作流程 (Workflow)

1. 从随机城镇开始。

2. 提议迁移到另一个随机选定的城镇。

3. 根据下列规则接受或拒绝：

   • 新城镇选民更多：总是接受
   • 新城镇选民更少：以概率接受

   $$\text{接受概率}=\min \left(1,\frac{\text{新城镇人数}}{\text{当前城镇人数}}\right)$$

4. 重复以上步骤大量次数（如 10000 次）。

5. 记录在各城镇停留的比例。

   举例说明：

   • 当前城镇人数是1000，新城镇人数是500：
     • 接受概率 = min(1, 500 / 1000) = 0.5
     • 随机抽U，如果U ≤ 0.5，就去新城镇；否则不去。
   • 当前城镇人数是300，新城镇人数是500：
     • 接受概率 = min(1, 500 / 300) = 1
     • 因为概率是1，直接去新城镇。

算法解释

• 类似于接受-拒绝方法。
• 利用马尔可夫链性质：每个新状态只依赖前一个状态（无记忆）。
• 目标是构建马尔可夫链 { X t , t = 0 , 1 , … } \{X_t, t=0,1,\ldots\} {Xt​,t=0,1,…}，使得极限分布为目标分布 $f(x)$。

算法泛化 (Generalisation)

• 从初始状态 $X_0$ 开始。
• 根据当前状态 $X_t$ 生成提议点 $Y$，服从建议分布 $q(y\mid x)$。
• 以接受概率 $\alpha (X_t,Y)$ 接受提议：

$$\alpha (x,y)=\min \left\{1,\frac{f(y)q(x\mid y)}{f(x)q(y\mid x)}\right\}$$

• 从均匀分布 $U\sim \text{Uniform}(0,1)$ 采样。
• 如果 $U\le \alpha$，则 $X_{t+1}=Y$，否则 $X_{t+1}=X_t$。

即使你的proposal distribution q(y|x) 不是完美的（比如不是精确的目标分布f(x)），通过Metropolis-Hastings算法这样设计的接受/拒绝机制，最终你仍然能够生成符合目标分布f(x)的样本。

换句话说：

• 你用一个“近似”的方式提出新点（proposal distribution），
• 通过计算接受概率 α 来修正这个近似，
• 这样长期下来，产生的样本序列会收敛到你想要的目标分布。

对称建议分布 (Symmetric Proposal)

• 如果建议密度是对称的，即

$$q(y\mid x)=q(x\mid y)$$

• 接受概率简化为

$$\alpha (x,y)=\min \left\{1,\frac{f(y)}{f(x)}\right\}$$

• 此时算法也称为随机游走采样器 (Random Walk Sampler)。
• 常用对称建议函数为高斯分布：

$$q(x)\sim \mathcal{N}(x_t,\sigma )$$

• $\sigma$ 控制状态空间的探索速度。

MCMC 的应用

贝叶斯统计中的后验采样

• 可以使用贝叶斯定理对模型和数据进行建模。

后验分布表示为：

$$p(\theta |\mathbf{D})=\frac{p(\mathbf{D}|\theta )p(\theta )}{p(\mathbf{D})}$$

其中：

• $p(\theta |\mathbf{D})$ 是参数 $\theta$ 在给定数据 $\mathbf{D}$ 后的后验概率。
• $p(\mathbf{D}|\theta )$ 是似然函数（likelihood），表示在参数为 $\theta$ 时，数据 $\mathbf{D}$ 出现的概率。
• $p(\theta )$ 是先验概率（prior），表示在获得数据前对参数的认识。
• $p(\mathbf{D})$ 是边缘似然（marginal likelihood）或证据，是归一化常数。

具体含义：

• Posterior：给定数据 $\mathbf{D}$ 后，参数 $\theta$ 发生的可能性。
• Likelihood ratio：数据 $\mathbf{D}$ 对参数 $\theta$ 的支持程度。
• Prior：未观察数据前参数 $\theta$ 的概率分布。

计算边缘似然通常是一个难以直接求解的积分：

$$p(\mathbf{D})=\int p(\mathbf{D}|\theta )p(\theta )d\theta$$

回归模型参数估计举例

考虑以下线性回归模型：

$$Y={\beta }_0+{\beta }_1{X}_1+\cdots +{\beta }_p{X}_p+\epsilon ,\epsilon \sim \mathcal{N}(0,{\sigma }^2)$$

参数集合为：

θ = { β 0 , β 1 , … , β p , σ 2 } \theta = \{\beta_0, \beta_1, \ldots, \beta_p, \sigma^2\} θ={β0​,β1​,…,βp​,σ2}

观察数据为：

$$\mathbf{D}=\{(Y_i,X_{i1},\dots ,X_{ip}),i=1,\dots ,n\}$$

似然函数基于正态分布，形式为 $p(\mathbf{D}|\theta )$，但后验分布：

$$p(\theta |\mathbf{D})\propto p(\mathbf{D}|\theta )p(\theta )$$

通常无法解析计算。此时，MCMC用于估计此后验分布，获取参数的分布特性。

计算难积分的间接方法

后验分布的归一化常数：

$$p(\mathbf{D})=\int p(\mathbf{D}|\theta )p(\theta )d\theta$$

难以直接计算积分，尤其在高维参数空间。

MCMC避免显式计算该积分，而是采样后验分布，间接获得参数分布信息。

抛硬币例子

假设我们观察多次抛硬币结果，目标是估计硬币正面概率：

$$\theta =P(\text{Head})$$

• 频率学派认为 $\theta$ 是一个固定但未知的常数；
• 贝叶斯派则视 $\theta$ 为随机变量，赋予先验分布，基于数据推断后验分布 $p(\theta |\text{Data})$。

MCMC用于从该后验分布采样，帮助估计硬币是否有偏。

采样相关性与Burn-in期

• MCMC开始阶段生成的样本，在算法收敛到真实分布之前，这段时间称为“burn-in period”（预热期）。

  • 这段样本应该被丢弃，不纳入最终分析。

• MCMC生成的样本是相关的，因为它们来源于一个马尔可夫链。

  • 以前，许多实践者主张对样本进行稀疏采样（thinning），即每隔第$k^{th}$个样本保留一个。
  • 这种做法历史上出于几个原因：
    • 减少样本间的相关性，更容易计算标准误差。
    • 减少存储马尔可夫链所需的空间。

总结

• MCMC为复杂概率分布提供有效采样方案，尤其适用于无法直接采样的后验分布。
• 广泛应用于贝叶斯统计参数估计、复杂模型推断、概率积分估计、事件概率推断等领域。
• 理解采样过程中的相关性及Burn-in期对结果质量至关重要。