【学习笔记】Matlab和python双语言的学习(蒙特卡洛法）

---

前言

通过模型算法，熟练对Matlab和python的应用。
学习视频链接：
https://www.bilibili.com/video/BV1EK41187QF?p=22&vd_source=67471d3a1b4f517b7a7964093e62f7e6

一、蒙特卡洛

• 蒙特卡洛严格意义上讲并不是一种算法模型，它是一种概率统计的思想方法，没有通用的代码，每个问题对应的代码都是不同的。
• 蒙特卡罗法又称统计模拟法，是一种随机模拟方法，以概率和统计理论方法为基础的一种计算方法，是使用随机数(或更常见伪随机数)来解决很多计算问题的方法。将所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这一方法的概率统计特征，故借用赌城蒙特卡罗命名。
• 只要求解的问题与概率模型有关联，就可以用蒙特卡洛法。

二、经典示例：计算圆周率 π

• 一个经典的蒙特卡洛方法的示例是计算圆周率 π。通过在单位正方形内随机生成点，估计圆的面积来逼近 π 的值。
• 步骤：
  1.在正方形[0,2]x[0,2] 内随机生成大量点。
  2.计算这些点中有多少点落在单位圆(圆心为 (1,1)，半径为 1)内。
  3.单位圆的面积是 π，正方形的面积是 4，因此可以通过以下公式估计 π:
  $$\pi \approx \frac{\text{圆内点的数量}}{\text{总点数}}\times 4$$

1.代码实现----Matlab

clear;clc
% 参数初始化：投放10000个点，圆的半径为1，圆心坐标(1，1)
% 初始时还未投放点，有0个点在圆内
p = 100000; r = 1; x0 = 1; y0 = 1; n = 0;
hold on
for i = 1:p
    px = rand*2;
    py = rand*2;
    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;

运行结果图：

2.代码实现----python

import matplotlib.pyplot as plt
import random
import numpy as np
# 计算圆周率
# 参数初始化：投放10000个点,圆的半径为1,圆心坐标(1，1)
# 初始时还未投放点,有0个点在圆内
p = 100000
r = 1
x0 = 1
y0 = 1
n = 0
# 设置绘图窗口
plt.figure()
plt.title('Monte Carlo Simulation for Esitimation Pi')
plt.xlabel('x')
plt.ylabel('y')
for i in range(p):
    px = np.random.rand() * 2
    py = np.random.rand() * 2
    if (px - x0) ** 2 + (py - y0) ** 2 < r ** 2:
        plt.plot(px,py,'.',color = 'b')
        n = n + 1
    else:
        plt.plot(px,py,'.',color = 'r')
plt.axis ('equal')
plt.show()

s = n / p * 4
pi0 = s

运行结果图：

三、示例2：三门问题

• 你参加一档电视节目，节目组提供了A、B、C三扇门，主持人告诉你，其中一扇门后边有辆汽车，其他两扇门后面是一头山羊，你可以选择一扇门打开获得门后的东西。
  假如你选择了B门，这时，主持人打开了C门，让你看到C门后是只山羊，然后问你要不要改选A门?(你想要汽车)

这个问题在概率论中是一个条件概率的问题，类似于在排除一个错误选项后，选择到正确选项的概率是多少。求解的问题与概率模型有关，我们使用蒙特卡洛法。

1.代码实现----Matlab

（1）有条件概率

clear;clc
% randi([a,b],m,n)  函数可在指定区间[a,b]内随机取出大小为m*n的整数矩阵
randi([1,5],5,8) % 在区间[1,5]内随机取出大小为5*8的整数矩阵

% 字符串的连接方式：(1)['字符串1','字符串2']  (2)strcat('字符串1','字符串2')

n = 100000;
a = 0;  % 不改变注意能赢的次数
b = 0;  % 改变注意能赢的次数
for i = 1:n
    x = randi([1,3]);    % 有车的那扇门
    y = randi([1,3]);   % 我们选择的那扇门
    if x == y
        a = a + 1; b = b + 0;    % 不改变主意能赢
    else
        a = a + 0; b = b + 1;    % 改变主意能赢
    end
end
disp(['使用蒙特卡洛方法得到不改变主意能赢的概率是', num2str(a/n)]);
disp(['使用蒙特卡洛方法得到改变主意才能赢的概率是', num2str(b/n)]);

运行结果：


（2）无条件概率（考虑没有获奖的情况）

clear;clc
% 考虑失败的情况
n = 100000;
a = 0;  % 不改变注意能赢的次数
b = 0;  % 改变注意能赢的次数
c = 0;  % 没有获奖的次数
for i = 1:n
    x = randi([1,3]);    % 有车的那扇门
    y = randi([1,3]);   % 我们选择的那扇门
    change = randi([0,1]);   % 改变主意是1，不改变主意是0
    if x == y   % 不改变才能赢
        if change == 0   % 不改变主意是0
            a = a + 1; b = b + 0; c = c + 0;  
        else
            a = a + 0; b = b + 0; c = c + 1;   
        end
    else        % 改变才能赢
        if change == 0   % 不改变主意是0
            a = a + 0; b = b + 0; c = c + 1;    
        else            % 改变主意是1
            a = a + 0; b = b + 1; c = c + 0;   
        end
    end
end
disp(['使用蒙特卡洛方法得到不改变主意能赢的概率是', num2str(a/n)]);
disp(['使用蒙特卡洛方法得到改变主意才能赢的概率是', num2str(b/n)]);
disp(['使用蒙特卡洛方法不获奖的概率是', num2str(c/n)]);

运行结果：

2.代码实现----python

（1）有条件概率

import numpy as np

n = 10000
a = 0   # 不改变主意能获胜的次数
b = 0   # 改变主意能获胜的次数

for i in range(n):
    x = np.random.randint(1,4)      #随机生成1-3的整数
    y = np.random.randint(1,4)
    if x == y:
        a += 1
    else:
        b += 1

print("使用蒙特卡洛方法得到不改变主意能赢的概率是", a/n)
print("使用蒙特卡洛方法得到改变主意能赢的概率是", b/n)

运行结果：

（2）无条件概率（考虑没有获奖的情况）

import numpy as np

n = 10000
a = 0   # 不改变主意能获胜的次数
b = 0   # 改变主意能获胜的次数
c = 0   # 不获胜的次数
for i in range(n):
    x = np.random.randint(1,4)      #随机生成1-3的整数
    y = np.random.randint(1,4)
    change = np.random.randint(0,2)  # % 改变主意是1，不改变主意是0
    if x == y:
        if change == 0:
            a += 1
        else:
            c += 1
    else:
        if change == 0:
            c += 1
        else:
            b += 1

print("使用蒙特卡洛方法得到不改变主意能赢的概率是", a/n)
print("使用蒙特卡洛方法得到改变主意能赢的概率是", b/n)
print("使用蒙特卡洛方法不能赢的概率是", c/n)

运行结果：

• 蒙特卡洛适用于概率统计的问题，但得出的结果只能用作估计，并不是一个精确的数值。代码中每次生成的随机数不同也会造成结果上的差异。使用蒙特卡洛法，当测试的数量足够大时，得出的结论才会更精确。

总结

本文介绍了使用蒙特卡洛法解决与概率统计有关联的问题。分析了两个示例，并分别使用Matlab和python进行代码编写。