蒙特卡罗介绍以及MATLAB算法实现 泡泡马特模拟抽盒

一、算法介绍

1. 背景

    蒙特· 卡罗方法(Monte Carlo method) ,也称统计模拟方法,是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明, 而被提出的一类以概率统计理论为指导的重要数值计算方法，其名字的灵感来源于摩纳哥的赌城“蒙特卡洛”。

2. 基本思想

    大家从小应该就知道，圆的面积公式: $\pi R^2$。但是这个 $\pi$ 得到的3.1415926…(位数太多我也记不住了qaq)是怎么得到的呢？接下来我将带大家进入蒙特卡洛的奇妙世界。

图1

    蒙特卡洛算法其实本质上很简单，Henry用一种更加简洁、直观的方式给大家介绍一下。
    假设在一个100×100（单位：$c{m}^2$）的空间中，各放置一个15×15的正方形盒子和半径为15的圆形盒子（此处不考虑盒子的高，因为对算法影响不大，可以假设高无线大，即可以容纳无限个乒乓球）。这时候，聪明的小明开始往这个空间里随机扔乒乓球，无意识地扔了100个、1000个、10000个…当数量足够大时，圆形盒子里乒乓球的个数和正方形盒子里的乒乓球个数的比值会趋向于一个值，而这个值，就是那个神奇的$\pi$。

3. 算法实现

1. 非动态算法演示

clc,clear,close all

square_side_length=30;
circle_D=square_side_length;
x1=10:10+square_side_length;
figure (1)
plot(x1,40*ones(square_side_length+1,1),'-k','LineWidth',2)
hold on
plot(10*ones(square_side_length+1,1),40:40+square_side_length,'-k','LineWidth',2)
plot(40*ones(square_side_length+1,1),40:40+square_side_length,'-k','LineWidth',2)
plot(x1,70*ones(square_side_length+1,1),'-k','LineWidth',2)
x2=50:110;
f_circle1=55+sqrt((square_side_length)^2-(x2-80).^2);
f_circle2=55-sqrt((square_side_length)^2-(x2-80).^2);
plot(x2,f_circle1,'-k','LineWidth',2)
plot(x2,f_circle2,'-k','LineWidth',2)
title('蒙特卡罗算法——小球随机下落(Henry.liu)','FontSize',14)
xlabel('X轴 单位：cm')
ylabel('Y轴 单位：cm')
axis([0 120 0 120])
set(gca,'Linewidth',2)
text(12,55,'落在正方形框里的小球','FontSize',6)
text(80,55,'落在圆框里的小球','FontSize',8)
Test_Number=10000;
count_square=0;
count_circle=0;
for i=1:Test_Number
    dot_x=120*rand;
    dot_y=120*rand;
    if 10<dot_x && dot_x<40 && 40<dot_y && dot_y<70
        count_square=count_square+1;
        plot(dot_x,dot_y,'r*','MarkerSize',1)
    elseif (dot_x-80)^2+(dot_y-55)^2<30^2
        count_circle=count_circle+1;
        plot(dot_x,dot_y,'bo','MarkerSize',2)
    else
        plot(dot_x,dot_y,'k*','MarkerSize',2)
    end
%     pause(0.001)
%     milliseconds
end

 pi_Estimate=count_circle/count_square;
 if abs(pi_Estimate-pi)/pi<=0.05
     disp(['圆和矩形面积比为Π,掉落比为',num2str(pi_Estimate),'，误差在5%以内，该实验证明了落点概率与面积成正比。'])
 else
     disp(['圆和矩形面积比为Π,掉落比为',num2str(pi_Estimate),'，误差在5%以外，该实验无法证明落点概率与面积成正比。'])
 end

2. 动态算法演示

clc,clear,close all

square_side_length=30;
circle_D=square_side_length;
x1=10:10+square_side_length;
figure (1)
plot(x1,40*ones(square_side_length+1,1),'-k','LineWidth',2)
hold on
plot(10*ones(square_side_length+1,1),40:40+square_side_length,'-k','LineWidth',2)
plot(40*ones(square_side_length+1,1),40:40+square_side_length,'-k','LineWidth',2)
plot(x1,70*ones(square_side_length+1,1),'-k','LineWidth',2)
x2=50:110;
f_circle1=55+sqrt((square_side_length)^2-(x2-80).^2);
f_circle2=55-sqrt((square_side_length)^2-(x2-80).^2);
plot(x2,f_circle1,'-k','LineWidth',2)
plot(x2,f_circle2,'-k','LineWidth',2)
title('蒙特卡罗算法——小球随机下落(Henry.liu)','FontSize',14)
xlabel('X轴 单位：cm')
ylabel('Y轴 单位：cm')
axis([0 120 0 120])
set(gca,'Linewidth',2)
text(12,55,'落在正方形框里的小球','FontSize',6)
text(80,55,'落在圆框里的小球','FontSize',8)
Test_Number=500;
count_square=0;
count_circle=0;
for i=1:Test_Number
    dot_x=120*rand;
    dot_y=120*rand;
    if 10<dot_x && dot_x<40 && 40<dot_y && dot_y<70
        count_square=count_square+1;
        plot(dot_x,dot_y,'r*','MarkerSize',1)
    elseif (dot_x-80)^2+(dot_y-55)^2<30^2
        count_circle=count_circle+1;
        plot(dot_x,dot_y,'bo','MarkerSize',2)
    else
        plot(dot_x,dot_y,'k*','MarkerSize',2)
    end
    pause(0.001)
    milliseconds
end

 pi_Estimate=count_circle/count_square;
 if abs(pi_Estimate-pi)/pi<=0.05
     disp(['圆和矩形面积比为Π,掉落比为',num2str(pi_Estimate),'，误差在5%以内，该实验证明了落点概率与面积成正比。'])
 else
     disp(['圆和矩形面积比为Π,掉落比为',num2str(pi_Estimate),'，误差在5%以外，该实验无法证明落点概率与面积成正比。'])
 end

（其实和非动态的区别就是中间加了一个pause()函数QAQ）
Henry批注：由于蒙特卡罗算法基于大量的尝试，所以大家如果想要更精确的结果，可以将代码中的 Test_Number的数值调大到尽可能大，只是代码运行时长可能会长一点，如果有需要也可以将plot代码删去，让程序只进行对乒乓球是否进盒子的判断，这样代码运行速度会更快。

二、具体实例

案例1：抽牌

    小树和小涵有一天晚上闲来无事，手里正好有一副牌，因此玩了一个游戏：小涵从这副牌中如果能连续抽中两张一样花色的牌，小树就做一个俯卧撑；若没抽中，自然小涵做一个俯卧撑。小涵心想，这不轻轻松松，于是便答应了下来。结果一晚上，小涵做了100多个俯卧撑，手都麻了，小树才做了二三十个。气不过的小涵连忙打开matlab，写下了下面这段代码，想试试到底自己这样的概率是多少：

 % 求一副卡牌（54张）中，一次最多摸两张，一下摸到两张花色相同牌的概率。
clc,clear,close all

Test_Number=1000000;
Possible_Test=zeros(1,Test_Number);
for i=1:Test_Number
   Random_Card=randperm(54);
   Full_Hand_red_peach=2;
   Top_Five=Random_Card(1,1:Full_Hand_red_peach);   
   if (sum(Top_Five<=13) ==2) || (sum(Top_Five>=14) ==2 &&  sum(Top_Five<=26) ==2) || (sum(Top_Five>=27) ==2 &&  sum(Top_Five<=39) ==2) || (sum(Top_Five>=40) ==2 &&  sum(Top_Five<=52) ==2)
       Possible_Test(i)=1;
   end
end
p=sum(Possible_Test)/Test_Number;
disp(['在一副牌中，一次最多摸两张，一下摸到两张',num2str(Full_Hand_red_peach),'张花色相同牌的概率是',num2str(p)])

    代码运算下来结果发现，摸到花色相同牌的概率仅有21%左右，小涵当场哭晕过去。

在一副牌中，一次最多摸两张，一下摸到两张2张花色相同牌的概率是0.21815

案例2：泡泡马特模拟抽盒

   有一天，小姜同学想给女朋友买59￥一盒的泡泡马特，他问小涵同学如果一次性抽中所有款大概要多少钱。小涵一想，兄弟别冲动啊！于是立马打开matlab进行了蒙特卡洛测试，测试代码如下所示：

clc,clear,close all

%% 抽齐所有款式所花的次数测试
Test_Number=10000;
type=[1*ones(1,13),2*ones(1,13),3*ones(1,13),4*ones(1,13),5*ones(1,13),6*ones(1,13),7*ones(1,13),8*ones(1,13),9*ones(1,13),10*ones(1,13),11*ones(1,13),12];

sum1=zeros(1,Test_Number);
for y=1:10000
    count1=0;
    total=zeros(1,12);
    i=1; 

    while 1
        x=randperm(144,1);
        m=type(x);
        if ~ismember(m,total)
            count1=count1+1;
            total(count1)=m;
            if count1==12
                break
            end
        end

        i = i + 1;
    end
    sum1(y)=i; % 记录每次蒙特卡洛试验抽中所有类型的泡泡马特所需要的次数
end
average_time=sum(sum1)/Test_Number
max(sum1)
min(sum1)
cost = average_time * 59

   经过计算，发现基本上要8800元左右才能将所有类型的泡泡马特集齐，究其原因还是因为泡泡马特的隐藏款太难抽了。小姜听到这个价格后，赶忙打消了这个念头。

average_time =

  149.1299


ans =

        1432


ans =

    14


cost =

   8.7987e+03

三、 总结

   蒙特卡洛算法，本质上是一种概率论的思想，但我们能用更直观的方式去理解他，概括下来就是一种能进行多次模拟的一种方法，不管是模拟上面的小球下落，还是后来的抽牌和抽盲盒，经过足够多次数的抽取，总能依概率收敛到某个数，这个数就是我们想得到的那个值。

ps: 由于中途经历了考研，所以本该是大概去年这个时候发的blog，一直拖到现在才发出来。后续会在数学建模专栏继续更新，也会结合本专业机器人方向开设新的内容，请大家多多关注！