原理

蒙特卡洛方法（Monte Carlo method）是一类用于解决各种问题的随机算法，广泛应用于数学、物理、工程、金融等领域。该方法以摩纳哥著名赌城蒙特卡洛命名，因为其核心思想依赖于随机抽样或统计抽样，类似于赌博中的随机性。

明确指令-拆解问题-分析生成

什么是蒙特卡洛法？

结论： 蒙特卡洛法是一种基于概率统计理论的数值计算方法，通过大量随机样本模拟来求解问题。

详细展开：

  1. 随机抽样： 蒙特卡洛方法利用随机数生成器产生大量的随机样本点，然后根据这些样本点来估计所研究的问题的解。
  2. 数值积分： 在无法解析地求出定积分时，可以通过蒙特卡洛方法来进行数值积分，即用随机点落在被积函数下的平均值乘以区域大小来近似积分值。
  3. 优化与仿真： 方法也常用于复杂系统的仿真和优化，如金融风险评估、粒子物理实验设计等。
  4. 误差估计： 随着样本数量的增加，蒙特卡洛方法给出的结果会逐渐收敛到真实值，且可以估算结果的不确定性或误差范围。

蒙特卡洛法的应用场景：

  • 物理科学： 模拟粒子传输、热力学性质等。
  • 金融工程： 计算期权定价、投资组合风险管理等。
  • 计算机图形学： 实现光线追踪，提高图像渲染的真实度。
  • 运筹学： 解决复杂的最优化问题，比如旅行商问题。

注意事项：

  • 计算效率： 对于高维问题，可能需要非常大的样本量才能获得足够精确的结果，这可能导致计算成本很高。
  • 随机数质量： 使用高质量的伪随机数生成器对于确保蒙特卡洛模拟的有效性至关重要。
  • 方差缩减技术： 为了提高精度并减少所需的样本数，可以采用诸如重要性抽样、控制变量等方差缩减技术。

代码

求π的近似值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
clc;clear
% 参数初始化：投放10000个点，圆半径为1，圆心坐标(1,1)
% 初始时还未投放点，有0个点在圆内
p = 10000;  r = 1; x0 = 1;  y0 = 1;  n = 0; 
hold on     % 保持绘图窗口，多次绘图
for i = 1:p     % 对于要投放的总共p个点
    % rand函数产生在(0, 1)之间的随机数；rand函数还有其他多种形式，可自行百度
    px = rand*2;    % 随机生成该点的横坐标
    py = rand*2;    % 随机生成该点的纵坐标
    % 所以，
    % 若该点在圆内，则颜色设为蓝色，变量n加一；在圆外则设为红色
    if (px-1)^2 + (py-1)^2 < 1      % 横纵坐标的平方和小于半径，则在圆内
        plot(px,py,'.','Color',"b");
        n = n+1;
    else
        plot(px,py,'.','Color',"r");
    end
end
axis equal      % 绘图时横纵坐标单位长度相同，便于观察圆
s = (n/p)*4;
pi0 = s

% 注意：matlab本身有圆周率值，在计算时直接调用pi即可
% a = 2*pi

1
pi0 = 3.1420

三门问题

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
%%  蒙特卡罗用于模拟三门问题
clear;clc
%% （1）预备知识
% randi([a,b],m,n)函数可在指定区间[a,b]内随机取出大小为m*n的整数矩阵
randi([1,5],5,8) %在区间[1,5]内随机取出大小为5*8的整数矩阵
%      2     5     4     5     3     1     4     2
%      3     3     1     5     4     2     1     2
%      4     1     3     3     2     2     5     1
%      5     3     3     4     4     5     4     4
%      4     2     3     4     2     4     2     4
randi([1,5])   %在区间[1,5]内随机取出1个整数
%     3

% 字符串的连接方式：(1)['字符串1','字符串2'] (2)strcat('字符串1','字符串2') 
['数模加油站','666']
strcat('数模加油站','666')

% num2str函数：将数值转换为字符串
mystr = num2str(1224)  % 注意观察工作区的mystr这个变量的值
disp([num2str(1224),'祝大家平安夜平平安安'])  % disp函数是输出函数

%% （2）代码部分（在成功的条件下的概率）
n = 100000;  % n代表蒙特卡罗模拟重复次数
a = 0;  % a表示不改变主意时能赢得汽车的次数
b = 0;  % b表示改变主意时能赢得汽车的次数
for i= 1 : n  % 开始模拟n次
    x = randi([1,3]);  % 随机生成一个1-3之间的整数x表示汽车出现在第x扇门后
    y = randi([1,3]);  % 随机生成一个1-3之间的整数y表示自己选的门
    % 下面分为两种情况讨论：x=y和x~=y
    if x == y   % 如果x和y相同，那么我们只有不改变主意时才能赢
        a = a + 1;     b = b + 0;
    else  % x ~= y ，如果x和y不同，那么我们只有改变主意时才能赢
        a = a + 0;     b = b +1;
    end
end
disp(['蒙特卡罗方法得到的不改变主意时的获奖概率为：', num2str(a/n)]);
disp(['蒙特卡罗方法得到的改变主意时的获奖概率为：', num2str(b/n)]);


%% （3）考虑失败情况的代码(无条件概率)
n = 100000;  % n代表蒙特卡罗模拟重复次数
a = 0;  % a表示不改变主意时能赢得汽车的次数
b = 0;  % b表示改变主意时能赢得汽车的次数
c = 0;  % c表示没有获奖的次数
for i= 1 : n  % 开始模拟n次
    x = randi([1,3]);  % 随机生成一个1-3之间的整数x表示汽车出现在第x扇门后
    y = randi([1,3]);  % 随机生成一个1-3之间的整数y表示自己选的门
    change = randi([0, 1]); % change =0  不改变主意，change = 1 改变主意
    % 下面分为两种情况讨论：x=y和x~=y
    if x == y   % 如果x和y相同，那么我们只有不改变主意时才能赢
        if change == 0  % 不改变主意
        	a = a + 1; 
        else  % 改变了主意
            c= c+1;
        end
    else  % x ~= y ，如果x和y不同，那么我们只有改变主意时才能赢
         if change == 0  % 不改变主意
        	c = c + 1; 
        else  % 改变了主意
            b= b + 1;
         end
    end
end
disp(['蒙特卡罗方法得到的不改变主意时的获奖概率为：', num2str(a/n)]);
disp(['蒙特卡罗方法得到的改变主意时的获奖概率为：', num2str(b/n)]);
disp(['蒙特卡罗方法得到的没有获奖的概率为：', num2str(c/n)]);

结果：

1
2
3
4
5
蒙特卡罗方法得到的不改变主意时的获奖概率为：0.33231
蒙特卡罗方法得到的改变主意时的获奖概率为：0.66769
蒙特卡罗方法得到的不改变主意时的获奖概率为：0.1672
蒙特卡罗方法得到的改变主意时的获奖概率为：0.33826
蒙特卡罗方法得到的没有获奖的概率为：0.49454