蒙特卡罗方法（Monte Carlo Method）是一大类随机算法的总称，是一种数值计算技术，其思想是其思想是通过随机样本来估计真实值：从问题的概率分布中生成大量随机样本，对每个样本进行模拟和计算，然后对所有样本的结果进行统计分析，以估计复杂系统的行为或特性。它广泛应用于积分计算、随机过程模拟、金融风险分析、物理系统模拟和优化问题等领域，依赖于大数定律和中心极限定理来确保结果的准确性和可靠性。

一、朴素蒙特卡罗（Crude Monte Carlo）方法

蒙特卡罗方法的核心思想是通过随机采样来模拟和估计复杂系统的行为。随机采样，简单粗暴、易于实现，是其最基本和直观的实现方式，因为它直接利用了概率和统计的基本原理。基于随机采样的蒙特卡罗方法称为朴素蒙特卡罗（Crude Monte Carlo）方法，以下是朴素蒙塔卡罗的一些应用场景：

1.1 蒙特卡罗方法近似计算$\pi$值

正方形面积：4
圆面积：pi
在正方形内随机投点n次，点落在圆内的概率：$p=\frac{\pi }{4}$，设落在圆内的点的样本数为m，则
$$p=\frac{m}{n}=\frac{\pi }{4}$$ $$\pi =\frac{4m}{n}$$

import numpy as np
import matplotlib.pyplot as plt

# 参数定义
n = 10000  # 样本数量
m = 0  # 落在圆中的样本数量

# 生成随机样本点
x = np.random.uniform(low=-1, high=1, size=n)
y = np.random.uniform(low=-1, high=1, size=n)

# 判断哪些点落在圆内
inside_circle = x**2 + y**2 <= 1

# 计算落在圆内的点的数量
m = np.sum(inside_circle)

# 近似计算 pi 的值
pi_approx = 4 * m / n
print("pi =", pi_approx)

# 可视化
plt.figure(figsize=(5, 5))
plt.scatter(x[inside_circle], y[inside_circle], color='blue', s=1, label='Inside Circle')
plt.scatter(x[~inside_circle], y[~inside_circle], color='red', s=1, label='Outside Circle')
plt.axhline(0, color='black', linewidth=0.5)
plt.axvline(0, color='black', linewidth=0.5)
plt.gca().set_aspect('equal', adjustable='box')
plt.title(f'Monte Carlo Approximation of Pi\nEstimated Pi = {pi_approx}')
plt.legend()
plt.show()

pi= 3.13912

\

1.2 蒙特卡罗方法近似定积分

一元函数定积分是指对一个变量的函数在某个区间上的积分。具体来说，给定一个函数 $f(x)$ 和一个区间 $[a,b]$，其定积分表示为：${\int }_a^{b}f(x)dx$。其几何意义是：函数$f(x)$在区间$[a,b]$上与x轴之间的面积。

1.2.1 蒙特卡罗方法近似定积分步骤

蒙特卡罗方法通过随机采样来近似计算一元函数定积分$I={\int }_a^{b}f(x)dx$。其基本步骤如下：

1. 在区间$[a,b]$做随机采样，得到n个样本点： { x i } i = 1 n \{x_i\}_{i=1}^n {xi​}i=1n​
2. 对每个样本点$x_i$，计算函数值$f(x_i)$
3. 计算这些函数值得平均值$\overline{f}=\frac{1}{N}{\sum }_{i=1}^{N}f(x_i)$
4. 估计积分值：${\int }_a^{b}f(x)dx\approx (b-a)\cdot \overline{f}$
5. 确定估计值的方差，通过找到估计的方差来量化我们的准确性，${\sigma }^2=\frac{b-a}{N}\sum _{i=1}^N{f}^2(x_i)-[\frac{b-a}{N}\sum _{i=1}^{N}f(x_i){]}^2$

1.2.2 示例

$f(x)=\frac{1}{1+\sin x\cdot (\ln x{)}^2}$，计算$f$在区间$[0.8,3]$的定积分：$I={\int }_{0.8}^{3}f(x)dx$

import numpy as np
import matplotlib.pyplot as plt

# 定义被积函数
f = lambda x: 1 / (1 + np.sin(x) * (np.log(x) ** 2))

# 蒙特卡罗积分参数
n = 1000
a, b = 0.8, 3

# 生成随机样本
x = np.random.uniform(low=a, high=b, size=n)

# 计算蒙特卡罗积分
I = (b - a) * np.mean(f(x))
print("蒙特卡罗积分结果:", I)

# 可视化
x_vals = np.linspace(a, b, 1000)
y_vals = f(x_vals)

plt.figure(figsize=(10, 6))
plt.plot(x_vals, y_vals, label='f(x)', color='black', linewidth=2)
plt.plot(x_vals, [np.mean(f(x))] * len(x_vals), label=r'$\bar{f}$', color='blue', linestyle='--', linewidth=2)

# 绘制蒙特卡罗估计的矩形
rect_height = np.mean(f(x))
plt.fill_between(x_vals, 0, [rect_height] * len(x_vals), color='blue', alpha=0.3)

# 设置图例和标签
plt.title('Monte Carlo Integration', fontsize=16)
plt.xlabel('x', fontsize=14)
plt.ylabel('f(x)', fontsize=14)
plt.legend(fontsize=12)
plt.grid(True, linestyle='--', alpha=0.7)

# 显示图形
plt.show()

蒙特卡罗积分结果: 1.7533011607192013

实际上用一个矩形估计定积分的面积。https://www.scratchapixel.com/lessons/mathematics-physics-for-computer-graphics/monte-carlo-methods-in-practice/monte-carlo-integration.html的这个图很好的展示了这个过程。

1.2.3 蒙特卡罗估计定积分优缺点

优点：

• 简单易行：实现简单，适用于高维积分和复杂区域，是一种强大的数值积分工具，尤其在高维和复杂区域的积分计算中具有独特优势。
• 适用广泛：不受函数形式和维数的限制。

缺点：

• 收敛速度慢：收敛速度较慢，通常需要大量样本点才能得到较高精度。
• 随机性：结果具有随机性，每次计算可能略有不同。

1.3 蒙特卡罗方法近似计算期望值

期望定义：
离散随机变量：如果 $X$ 是一个离散随机变量，且其取值为 $x_1,x_2,\dots ,x_n$，对应的概率为 $P(X=x_i)=p_i$，则期望值 $E(X)$ 定义为：$E(X)={\sum }_{i=1}^n{x}_i\cdot p_i$。
连续随机变量：如果 $X$ 是一个连续随机变量，且其概率密度函数为$f(x)$，则期望值 $E(X)$定义为：$E(X)={\int }_{-\infty }^{\infty }x\cdot f(x)dx$

1.3.1 蒙特卡罗方法近似计算期望值步骤

蒙特卡罗方法通过随机采样来近似计算期望值。其基本步骤如下：

1. 确定随机变量$X$及其分布
2. 根据随机变量X的分布随机采样，生成N个样本$x_1,x_2,\dots ,x_N$
3. 计算样本的平均值$\overline{X}=\frac{1}{N}{\sum }_{i=1}^N{x}_i$
4. 估计期望，样本均值$\overline{X}$即为期望值的近似估计

例子
好的，让我们通过一个具体的例子来说明如何使用蒙特卡洛方法来近似计算均值。我们将计算一个简单函数在特定区间上的均值。
例子：计算函数 f(x) = x^2 在区间 [0, 1] 上的均值
步骤 1: 定义问题
我们想要计算函数 f(x) = x^2 在区间 [0, 1] 上的均值。理论上，这个均值可以通过积分来计算：
E[f(x)] = \int_0^1 x^2 , dx
这个积分的结果是 \frac{1}{3}。
步骤 2: 选择概率分布
在这个问题中，我们可以直接从区间 [0, 1] 上的均匀分布中采样，因为函数的定义域就是这个区间。
步骤 3: 生成随机样本
我们生成 n 个随机样本 x_1, x_2, …, x_n，每个样本都是从 [0, 1] 区间上均匀分布的随机数。
步骤 4: 计算函数值
对于每个随机样本 x_i，计算 f(x_i) = x_i^2。
步骤 5: 计算样本均值
计算这些函数值的样本均值：
\hat{E}[f(x)] = \frac{1}{n} \sum_{i=1}^{n} f(x_i) = \frac{1}{n} \sum_{i=1}^{n} x_i^2
步骤 6: 评估误差
随着 n 的增加，样本均值 \hat{E}[f(x)] 将越来越接近真实的均值 \frac{1}{3}。
实际计算
假设我们生成了 1000 个随机样本，以下是 Python 代码来实现这个过程：

import numpy as np

# 设置随机种子以确保结果可重复
np.random.seed(0)

# 生成1000个从[0, 1]均匀分布的随机样本
samples = np.random.uniform(0, 1, 1000)

# 计算函数值
function_values = samples**2

# 计算样本均值
sample_mean = np.mean(function_values)

print("近似的均值:", sample_mean)

运行这段代码，你将得到一个接近 \frac{1}{3} 的结果，即大约 0.333。
这个例子展示了如何使用蒙特卡洛方法来近似计算一个简单函数在特定区间上的均值。通过增加样本数量，你可以进一步提高近似的精度。

1.3.3 示例

$$f(x)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{x^2}{2}}$$

1.5 蒙特卡罗采用应用于图像处理

蒙特卡罗方法是一种统计学方法，广泛应用于数值积分、优化和其他计算问题。对于图像处理，蒙特卡罗方法可以用于图像采样、降噪、光照模拟等任务。下面是一个简单的示例，展示如何使用蒙特卡罗方法对图像进行采样。

import random
from PIL import Image, ImageDraw


def monte_carlo_sampling(image_path, num_samples):
    # 打开图像
    image = Image.open(image_path)
    width, height = image.size

    # 创建一个新的图像来绘制采样点
    sampled_image = Image.new("RGB", (width, height), (255, 255, 255))
    draw = ImageDraw.Draw(sampled_image)

    # 使用蒙特卡罗方法进行采样
    for _ in range(num_samples):
        x = random.randint(0, width - 1)
        y = random.randint(0, height - 1)
        color = image.getpixel((x, y))
        draw.point((x, y), fill=color)

    return sampled_image


# 示例使用
image_path = "./蒙特卡罗/R-C.jpg"
num_samples = 1000000  # 采样点的数量
sampled_image = monte_carlo_sampling(image_path, num_samples)

# 保存或显示采样后的图像
sampled_image.show()
sampled_image.save("./蒙特卡罗/sampled_image.jpg")

原图（左），采样（图右）

蒙特卡罗方法在图像处理中的应用场景非常广泛，主要利用其随机采样和统计特性来解决一些复杂的问题。以下是一些具体的应用场景：

• 图像降噪：蒙特卡罗方法可以用于图像降噪，通过多次采样和平均来减少噪声。例如，在光线追踪渲染中，蒙特卡罗方法可以用于模拟光线的传播路径，从而生成更逼真的图像。
• 光照模拟： 在计算机图形学中，蒙特卡罗方法常用于模拟光照和阴影效果。通过随机采样光线的路径，可以逼真地模拟光线在场景中的传播和反射，从而生成高质量的渲染图像。
• 图像重建： 蒙特卡罗采样可以用于从部分数据重建图像。例如，在医学成像中，可以通过随机采样的方式重建缺失的图像数据。
• 纹理生成： 蒙特卡罗方法可以用于生成随机纹理，通过随机采样和合成来创建自然的纹理效果，如云、烟雾、水波等。
• 图像压缩： 蒙特卡罗采样可以用于图像压缩，通过随机选择重要的像素点进行存储，从而减少数据量。
• 图像分割和特征提取： 蒙特卡罗方法可以用于图像分割和特征提取，通过随机采样像素点来识别和分割图像中的不同区域。
• 机器学习和计算机视觉： 在机器学习和计算机视觉中，蒙特卡罗方法可以用于数据增强和模型训练，通过随机采样生成多样化的训练数据，提高模型的泛化能力。

二、朴素蒙特卡罗（Crude Monte Carlo）的局限性

朴素蒙特卡罗（Crude Monte Carlo）方法是一种简单直接的随机采样方法，尽管它在许多应用中非常有用，但也存在一些显著的缺点和局限性。以下是朴素蒙特卡罗方法的一些主要缺点：

1. 低效率
   • 高方差：朴素蒙特卡罗方法通常具有较高的方差，这意味着需要大量的样本才能获得精确的估计。对于高维问题，所需的样本数量可能会非常大，从而导致计算成本高昂。
   • 慢收敛：由于高方差，朴素蒙特卡罗方法的收敛速度较慢，需要大量的计算资源和时间来达到所需的精度。
2. 维度灾难。高维问题：在高维空间中，朴素蒙特卡罗方法的效率显著下降。随着维度的增加，样本点在空间中的分布变得稀疏，导致估计的精度显著降低。这种现象被称为“维度灾难”。
3. 采样不均匀。不均匀采样：朴素蒙特卡罗方法在整个采样空间内均匀地生成样本，但在许多实际问题中，目标函数在某些区域的变化更为显著。均匀采样可能会导致在重要区域的样本不足，从而影响估计的精度。
4. 计算成本高
   大量样本：为了获得精确的估计，朴素蒙特卡罗方法通常需要生成大量的样本，这会导致计算成本和存储需求显著增加。
5. 缺乏自适应性。固定采样策略：朴素蒙特卡罗方法使用固定的采样策略，无法根据目标函数的特性动态调整采样分布。这意味着它不能有效地利用已有的样本信息来改进采样策略。
6. 不适用于复杂分布。复杂分布：对于具有复杂形状或多峰的目标分布，朴素蒙特卡罗方法可能无法有效地采样，导致估计结果不准确。
7. 难以处理约束条件。约束条件：在处理具有约束条件的问题时，朴素蒙特卡罗方法可能需要额外的步骤来确保生成的样本满足约束条件，这会进一步增加计算复杂度。
8. 缺乏方向性。无方向性：朴素蒙特卡罗方法没有利用任何关于目标函数的梯度或其他方向性信息，这使得它在优化问题中效率较低。

三、蒙特卡罗改进方法

尽管简单性和广泛适用性使其成为解决许多复杂问题的有力工具，但其在某些情况下效率较低。在实际应用中，可以结合其他改进方法，如重要性采样、马尔可夫链蒙特卡罗（MCMC）等，以提高采样效率和估计精度。

• 3.1 重要性采样（Importance Sampling）：通过改变采样分布，使得更多的样本来自对结果贡献较大的区域，从而减少方差，提高估计的精度。

• 3.2 自适应重要性采样（Adaptive Importance Sampling：动态调整采样分布，使其逐渐逼近最优分布，从而进一步提高采样效率。

• 3.3 拒绝采样（Rejection Sampling：通过在一个简单的分布上生成样本，并根据目标分布的概率密度函数进行筛选，保留符合条件的样本。

• 3.4 马尔可夫链蒙特卡罗（MCMC）方法
  通过构建一个马尔可夫链来生成样本，样本的分布逐渐逼近目标分布。常见的MCMC方法包括：

• Metropolis-Hastings算法：通过接受-拒绝机制生成样本。

• Gibbs采样：逐步更新每个变量的条件分布来生成样本。

• 3.5 拉丁超立方采样（Latin Hypercube Sampling：将每个维度的采样空间分成等概率的区间，然后在每个区间内随机采样，确保每个维度的样本均匀分布。

• 3.6 低差异序列（Quasi-Monte Carlo：使用低差异序列（如Sobol序列、Halton序列）来进行采样，减少方差，提高积分估计的精度。

• 3.7 多重重要性采样（Multiple Importance Sampling：结合多个重要性采样分布，通过加权平均来生成样本，进一步减少方差。

• 3.8 分层采样（Stratified Sampling：将采样空间划分为若干子区域，在每个子区域内独立采样，确保每个子区域的样本均匀分布，减少方差。

• 3.8 重要性重采样（Resampling：在粒子滤波等应用中，通过对已有样本进行重采样，保留高权重的样本，丢弃低权重的样本，提高样本质量。

• 3.10 变分推断（Variational Inference：通过优化一个简单的分布，使其逼近目标分布，从而生成高质量的样本。常用于贝叶斯推断和深度学习。

• 3.11 平行蒙特卡罗（Parallel Monte Carlo：利用并行计算资源，同时生成多个样本，提高采样效率。常用于高性能计算和大规模模拟。

• 3.12 多重链蒙特卡罗（Multiple Chain Monte Carlo：同时运行多个马尔可夫链，通过交换样本和信息，提高采样效率和收敛速度。

• 3.13 Sequential Monte Carlo（SMC）方法
  也称为粒子滤波，通过一系列的时间步长来更新样本分布，常用于动态系统的状态估计。

四、样本均值近似算法在随机规划问题中的应用

4.1 样本均值近似算法解决随机规划的思想

样本均值近似方法简称（Sample Average Approximation，SAA），是蒙特卡罗方法的一种具体应用，对一类具有随机分布特性的最优解进行了研究。其原理是先产生一个随机样点，再用样点的均值来逼近期望值函数。对所得到的样本均值逼近问题展开求解，这个过程被反复多次，直到估计的目标值与下界的差值低于某一阈值时，才会停止迭代，从而获得满意的解。

假设不确定需求值服从正态分布，有确定的分布函数，可以采用样本均值近似算法，对随机变量x进行n次随机抽取样本，生成独立同分布的样本，称为不同情景s。此时不确定性问题即转化为若干个确定性优化问题，不同情景下的样本均值即为所求解。

不同情景的作用在于将不确定性问题转化为若干个确定性优化问题。每个情景代表一种可能的需求情况，通过生成多个可能的情景（scenarios），在每个情景下进行确定性优化，从而在不确定性环境中找到一个稳健的解决方案。

4.2 理论基础

采用样本均值近似算法的核心思想是利用大数定律和中心极限定理，这些定理为我们提供了理论基础，使得在有确定分布的情况下，样本均值可以很好地近似总体的期望值。以下是详细解释：
大数定律
大数定律（Law of Large Numbers）是概率论中的一个基本定理，它表明，当样本数量 $n$ 足够大时，样本均值会收敛到总体的期望值。具体来说，如果 $X_1,X_2,\dots ,X_n$ 是从总体中抽取的独立同分布的随机变量，且每个随机变量的期望值为 $\mu$，则样本均值 $\overline{X}$ 会收敛到 $\mu$：
$$\overline{X}=\frac{1}{n}\sum _{i=1}^n{X}_i\to \mu \text{当}n\to \infty$$

这意味着，通过增加样本数量，我们可以使样本均值越来越接近总体的期望值。

中心极限定理
中心极限定理（Central Limit Theorem）进一步支持了样本均值近似算法的有效性。该定理表明，对于一组独立同分布的随机变量 $X_1,X_2,\dots ,X_n$，无论这些随机变量的原始分布是什么，当样本数量 $n$ 足够大时，样本均值的分布将近似于正态分布，即：$\overline{X}\sim N\left(\mu ,\frac{{\sigma }^2}{n}\right)$

其中，$\mu$是总体的期望值，${\sigma }^2$ 是总体的方差。这意味着，即使原始数据不是正态分布，样本均值的分布在大样本情况下也会接近正态分布。

为什么有确定的分布，就可以采用样本均值近似算法？

• 期望值的存在：有确定的分布意味着我们知道总体的概率分布函数 ( f(x) )，从而可以确定总体的期望值 ( \mu ) 和方差 ( \sigma^2 )。这为我们提供了一个明确的目标，即通过样本均值来估计这个期望值。
• 独立同分布样本：从确定的分布中抽取的样本是独立同分布的，这满足大数定律和中心极限定理的前提条件。
• 样本均值的收敛性：根据大数定律，当样本数量 $n$ 足够大时，样本均值 $\overline{X}$ 会收敛到总体的期望值 $\mu$。这使得样本均值成为一个可靠的估计量。
• 样本均值的分布特性：根据中心极限定理，当样本数量 $n$ 足够大时，样本均值的分布近似于正态分布。这使得我们可以使用正态分布的性质来进行进一步的统计推断和置信区间估计。

参考

• https://blog.csdn.net/robert_chen1988/article/details/90319257
• https://blog.csdn.net/wdl1992/article/details/108410404
• https://zhuanlan.zhihu.com/p/39628670
• https://zhuanlan.zhihu.com/p/300928400
• https://zhuanlan.zhihu.com/p/110855802