5.1.1 蒙特卡洛

什么是蒙特卡洛？

蒙特卡洛（Monte Carlo）是一种思想或一类方法，它通过随机采样来估计数学问题的数值解，本质上是通过大量采样来暴力解决问题，特别适用于复杂的积分、概率计算和优化问题等。

例子：如何利用蒙特卡洛方法估计圆周率 $\pi$ ？

问题描述：在一个边长为 2 的正方形中，内嵌一个半径为 1 的圆，我们希望估计圆的面积。
生成随机样本：在正方形内随机生成 $n$ 个点，记录落在圆内的点数 $m$
计算估计值：
- 圆的面积与正方形面积的比值是 $\frac{\pi r^2}{(2r)^2} = \frac{\pi}{4}$
- 通过 $\frac{m}{n}$ 估计这个比值， $\pi \approx 4 \times \frac{m}{n}$
增加样本数量：随着样本数量 $n$ 增加， $\pi$ 的估计值会越来越接近真实值

在未知 $\pi$ 的情况下，我们无法直接计算圆的面积，于是通过暴力采样和计算点数量的方式估计圆的面积，进而估计 $\pi$ 的值。

蒙特卡洛如何帮助我们处理复杂后验分布？

在求解后验分布的时候，我们经常要在分布上做一些计算，比如求期望、方差、或其他积分。但是如果碰上复杂分布、高维变量，就会让计算变得特别困难。这时候我们就可以尝试绕过解析解，通过大量采样的方法直接得到统计量。

例如，我们想对一个连续随机变量 $X$ 求期望，概率分布函数为：

p(x) = e^{-0.5 \cdot x^2} + \sin(x) + x^3

这是一个非标准概率密度函数，我们无法直接求积分，但是我们可以从 $p(x)$ 中采样 $x_1, x_2, x_3, ...x_i$ ，再利用大数定理（在本书中不做详细解释）近似期望：

\mathbb{E}[X] = \int x_i\cdot p(x_i)\,d(x)\approx \frac{1}{n} \sum_{i}^{n} x_i

同样，对离散连续变量 $X$ 求期望，虽然不涉及积分困难，但当变量的维度很高的时候 $(X \in \mathbb{R}^D )$ ，会出现维数灾难，求和计算难度随着维度增加而指数级增加。我们也可以通过从 $P(x)$ 采样来近似期望：

\mathbb{E}[X] = \sum_{i}^{n} x_i\cdot P(x_i)\approx \frac{1}{n} \sum_{i}^{n}x_i

一句话总结：遇到棘手的概率分布形式，避开难求的解析解，通过大量采样进行暴力计算估计。

注意：这里我们已知概率函数 $p(x)$ , 只是这个 $p(x)$ 可能比较复杂，直接积分比较难求。

蒙特卡洛的实际应用案例

接下来我们借着一个模拟股票走势的代码案例，展示蒙特卡洛的用法。

假定某个股票目前的股价是 $S_t$ , 预测时间 $\Delta t$ 之后的股价 $S_{t+1}$ ，金融工程中有公式：

S_{t+1} = S_t + \hat{\mu} S_t \Delta t + \sigma S_t \epsilon \sqrt{\Delta t}

假定初始股价 $S_0=10$ ，我们可以模拟一下，一年之后 (244 个交易日 ) 的股价分布。

首先，定义初始状态和模拟次数。

from scipy.stats import norm
import numpy as np
import matplotlib.pyplot as plt

s0 = 10. #刚买入的股价，初始状态
mu = 0.15
sigma = 0.2 
nSimu = 10000 # 模拟多少次
nStep = 244 # 每一次模拟走多少步, 1年只有244个交易日
dt = 1/nStep # 由于公式中的单位是年，而一年中有244个交易日，所以dt是理想化的每次股价变化间隔时间。

模拟10000次，得到一年后的股价分布。

price = np.empty((nSimu, nStep))
e = norm()
for i in range (nSimu):
    s = s0
    price[i, 0] = s
    for j in range(nStep):
        et = e.rvs() # 
        s = s + mu*s*dt + sigma*s*et*np.sqrt(dt) #股价的计算公式
        price[i, j] = s

画出一年后股价的分布。

plt.hist(price[:,-1], bins=30, density=True)
plt.ylabel('pdf')
plt.xlabel('Price')

最后更新于1天前