蒙特卡洛抽样方法详解

蒙特卡洛抽样方法详解

引用

CSDN

1.

https://m.blog.csdn.net/weixin_46713695/article/details/140424165

蒙特卡洛抽样方法是一种通过随机采样来进行数值计算的方法。它常用于求解复杂的数学问题，尤其是那些难以通过解析方法直接求解的问题。本文将详细介绍蒙特卡洛抽样方法的基本概念、具体步骤、应用场景以及代码实现。

认识该方法

蒙特卡洛抽样方法的核心思想是通过大量的随机样本来逼近一个确定的数值结果。例如，计算π的值时，可以通过在正方形内随机打点，根据落在1/4圆内的点的比例来估算π的值。抽样点越多，估算结果越准确。

详细介绍

基本步骤

1. 定义问题：确定需要估计的目标量。例如，积分值、概率值或某个随机过程的特性。
2. 构建模型：建立一个能够生成随机样本的数学模型。这个模型应能够代表问题的统计特性。
3. 生成随机样本：使用随机数生成器从模型中抽取大量随机样本。
4. 计算统计量：对随机样本进行计算，以得到所需的统计量，例如均值、方差等。
5. 逼近结果：根据大量样本的统计量来估计目标量。

应用

1. 积分计算：对于高维度的积分问题，蒙特卡洛方法可以通过随机采样来估计积分值。
2. 金融领域：用于期权定价、风险评估等。例如，蒙特卡洛模拟可以用于估算金融资产的未来价格分布。
3. 物理学和工程学：模拟粒子运动、热传导、辐射传输等复杂物理过程。
4. 优化问题：在求解复杂的优化问题时，可以使用蒙特卡洛方法来搜索最优解。

优点和缺点

优点：

• 能够处理高维度问题。
• 对模型的假设较少，灵活性强。
• 简单易实现，适用于各种复杂问题。

缺点：

• 收敛速度较慢，需要大量样本才能得到精确结果。
• 结果的精度依赖于样本数量，计算成本较高。

示例

假设我们需要计算单位圆内的面积，可以通过以下步骤使用蒙特卡洛抽样来估计：

1. 在正方形[0, 1] × [0, 1]内随机生成大量点。
2. 计算有多少点落在单位圆内（即点(x, y)满足x² + y² ≤ 1）。
3. 用落在圆内的点数与总点数的比值乘以正方形的面积（即1）来估计圆的面积。

import random

def monte_carlo_pi(num_samples):
    inside_circle = 0
    for _ in range(num_samples):
        x, y = random.uniform(0, 1), random.uniform(0, 1)
        if x**2 + y**2 <= 1:
            inside_circle += 1
    return (inside_circle / num_samples) * 4

# 估算圆周率
estimated_pi = monte_carlo_pi(1000000)
print(f"Estimated Pi: {estimated_pi}")

上述代码通过生成1000000个随机点来估算圆周率，结果会接近3.14159。

多维度随机抽样应用

代码逻辑解析

下面这段代码使用蒙特卡洛抽样方法来生成某些随机变量的样本。它首先定义了一个输入分布，然后从该分布中生成随机样本。下面逐步解析代码的逻辑：

1. 定义输入分布

input_dist_uncoupled = [
    ('v50', ot.TruncatedDistribution(v50_dist, 25., 50.)),
    ('air_den', ot.TruncatedDistribution(rou_dist, 0.8, 1.35)),
    ('ti', ot.Uniform(0.025, 0.2)),
    ('shear', ot.TruncatedDistribution(shear_dist, 0, 0.4)),
]

• 这里定义了四个输入变量的分布：
• v50：截断分布（范围在25到50之间）。
• air_den：截断分布（范围在0.8到1.35之间）。
• ti：均匀分布（范围在0.025到0.2之间）。
• shear：截断分布（范围在0到0.4之间）。
  这些分布使用了openturns库中的TruncatedDistribution和Uniform。

2. 创建组合分布

distribution_uncoupled = ot.ComposedDistribution([dist[1] for dist in input_dist_uncoupled])
distribution_uncoupled.setDescription([dist[0] for dist in input_dist_uncoupled])

• ot.ComposedDistribution用于创建一个组合分布，它包含了上述定义的所有分布。
• setDescription用于给分布中的每个变量命名，便于后续处理。

3. 生成蒙特卡洛样本

n_sample_mc = 100
ot.RandomGenerator.SetSeed(0)
sample_mc = np.array(distribution_uncoupled.getSample(n_sample_mc))
input_sample_monte_carlo = pd.DataFrame(data=sample_mc, columns=distribution_uncoupled.getDescription())

• n_sample_mc = 100：定义了生成样本的数量，这里设置为100。
• ot.RandomGenerator.SetSeed(0)：设置随机数生成器的种子，以确保结果可重复。
• distribution_uncoupled.getSample(n_sample_mc)：从组合分布中生成n_sample_mc个样本。
• pd.DataFrame(data=sample_mc, columns=distribution_uncoupled.getDescription())：将生成的样本转换为一个DataFrame，列名为之前定义的变量名。

代码整体逻辑

1. 定义输入变量的概率分布：确定每个输入变量的概率分布，包括截断范围和均匀分布。
2. 创建组合分布：将所有输入变量的分布组合在一起，形成一个多维的组合分布。
3. 生成随机样本：使用蒙特卡洛方法从组合分布中生成指定数量的随机样本。
4. 组织数据：将生成的样本数据整理为一个DataFrame，方便后续分析和处理。

应用场景

这段代码可以用于模拟不同的输入变量组合对某个系统或模型的影响。例如，在风能领域，可以模拟不同的风速（v50）、空气密度（air_den）、湍流强度（ti）和风切变（shear）对风力发电机性能的影响。通过蒙特卡洛方法生成的样本，可以进行不确定性分析、敏感性分析或风险评估等。

详述使用蒙特卡洛方法从组合分布中生成指定数量的随机样本的具体逻辑

1. 定义输入分布

首先，我们定义了每个输入变量的分布。这包括风速(v50)、空气密度(air_den)、湍流强度(ti)和风切变(shear)的概率分布：

input_dist_uncoupled = [
    ('v50', ot.TruncatedDistribution(v50_dist, 25., 50.)),
    ('air_den', ot.TruncatedDistribution(rou_dist, 0.8, 1.35)),
    ('ti', ot.Uniform(0.025, 0.2)),
    ('shear', ot.TruncatedDistribution(shear_dist, 0, 0.4)),
]

在这段代码中：

• v50_dist是v50的基础分布，被截断到[25, 50]范围内。
• rou_dist是air_den的基础分布，被截断到[0.8, 1.35]范围内。
• ot.Uniform(0.025, 0.2)定义了ti的均匀分布。
• shear_dist是shear的基础分布，被截断到[0, 0.4]范围内。

2. 创建组合分布

我们将这些输入变量的分布组合成一个多维的组合分布：

distribution_uncoupled = ot.ComposedDistribution([dist[1] for dist in input_dist_uncoupled])
distribution_uncoupled.setDescription([dist[0] for dist in input_dist_uncoupled])

在这段代码中：

• ot.ComposedDistribution([dist[1] for dist in input_dist_uncoupled])创建一个组合分布，包含了所有单变量分布。
• setDescription([dist[0] for dist in input_dist_uncoupled])为组合分布中的每个变量设置名称。

3. 生成蒙特卡洛样本

接下来，我们使用蒙特卡洛方法从组合分布中生成随机样本：

n_sample_mc = 100
ot.RandomGenerator.SetSeed(0)
sample_mc = np.array(distribution_uncoupled.getSample(n_sample_mc))
input_sample_monte_carlo = pd.DataFrame(data=sample_mc, columns=distribution_uncoupled.getDescription())

详细步骤如下：

1. 定义样本数量：

n_sample_mc = 100

这里，我们设置要生成的样本数量为100。

2. 设置随机数种子：

ot.RandomGenerator.SetSeed(0)

设置随机数生成器的种子为0，以确保生成的随机样本是可重复的。这对于调试和验证结果非常重要。

3. 生成随机样本：

sample_mc = np.array(distribution_uncoupled.getSample(n_sample_mc))

使用getSample(n_sample_mc)方法从组合分布中生成100个样本。这个方法会根据组合分布中的每个变量的定义，随机抽取100组变量值。生成的样本是一个二维数组，每行代表一个样本，每列代表一个变量。

4. 转换为DataFrame：

input_sample_monte_carlo = pd.DataFrame(data=sample_mc, columns=distribution_uncoupled.getDescription())

将生成的样本数据转换为一个Pandas DataFrame，并设置列名为之前定义的变量名称。这使得数据更易于管理和分析。

具体应用流程图

1. 定义各输入变量的概率分布：

• 使用ot.TruncatedDistribution和ot.Uniform定义各变量的分布。
• 设定截断范围和均匀分布范围。

2. 创建组合分布：

• 使用ot.ComposedDistribution将各变量的分布组合成一个多维分布。
• 使用setDescription为每个变量设置名称。

3. 生成蒙特卡洛样本：

• 设定要生成的样本数量。
• 设置随机数种子以确保结果可重复。
• 使用getSample方法从组合分布中生成随机样本。
• 将生成的样本数据转换为Pandas DataFrame，并设置列名。

通过上述流程，我们可以从定义的组合分布中生成指定数量的随机样本，并将其组织为易于分析的数据格式。这些样本可以用于进一步的统计分析、建模和仿真。