蒙特卡洛思想

算法思想

本质就是通过类似抽样调查的方法对于未知量进行判断结果，主要针对离散问题的仿真，有时可以实现连续性问题向离散型问题的转化，进而用蒙特卡洛进行估计

蒲丰投针

l =  0.520;    			 % 针的长度（任意给的）
a = 1.314;    			% 平行线的宽度(大于针的长度l即可)
n = 1000;    			% 做n次投针试验，n越大求出来的pi越准确
m = 0;    				% 记录针与平行线相交的次数
x = rand(1, n) * a / 2 ;   		% 在[0, a/2]内服从均匀分布随机产生n个数， x中每一个元素表示针的中点和最近的一条平行线的距离
phi = rand(1, n) * pi;    		% 在[0, pi]内服从均匀分布随机产生n个数，phi中的每一个元素表示针和最近的一条平行线的夹角
axis([0,pi, 0,a/2]);   box on;  % 画一个坐标轴的框架，x轴位于0-pi，y轴位于0-a/2， 并打开图形的边框
for i=1:n  						% 开始循环，依次看每根针是否和直线相交
    if x(i) <= l / 2 * sin(phi (i))     % 如果针和平行线相交
        m = m + 1;    					% 那么m就要加1
        plot(phi(i), x(i), 'r.') 		% 模仿书上的那个图，横坐标为phi，纵坐标为x , 用红色的小点进行标记
        hold on  						% 在原来的图形上继续绘制
    end
end
p = m / n;   							 % 针和平行线相交出现的频率
mypi = (2 * l) / (a * p); 				 % 我们根据公式计算得到的pi
disp(['蒙特卡罗方法得到pi为：', num2str(mypi)])

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from random import *
import warnings
warnings.filterwarnings("ignore")
L = 1
A = 2
num = 100000
res = 0
x = np.random.rand(num) * A / 2
phi = np.random.rand(num) * np.pi
fig=plt.figure()
ax=fig.add_subplot(111)
for i in range(0, num):
    if x[i] <= L / 2 * math.sin(phi[i]):
        res = res + 1
        ax.scatter(phi[i],x[i],c="r")
p = res / num
my_pi = (2 * L) / (A * p)
print(my_pi)
ax.set_xlim([0,math.pi])
ax.set_ylim([0,A/2])
plt.show()

三门问题

all=10000;
change=0;
notchange=0;
for i = 1:all
    r=randi([1,3]);
    tr=randi([1,3]);
    if r==tr
        notchange=notchange+1;
    else
        change=change+1;
    end
end
disp(["不改变：",notchange/all]);
disp(["改变：",change/all])

排队问题

arrive(1)=exprnd(10);
start(1)=arrive(1);
over(1)=normrnd(10,2)+start(1);
i=2;
timeSum=0;
while true
    arrive(i)=arrive(i-1)+exprnd(10);
    start(i)=max(arrive(i),over(i-1));
    if start(i)>480
        break
    end
    timeSum=timeSum+(start(i)-arrive(i));
    thisTime=normrnd(10,2);
    if thisTime<1
        over(i)=start(i)+1;
    else
        over(i)=start(i)+thisTime;
    end
    i=i+1;
end
disp(["一天共有"+num2str(i-1)+"人"])
disp(["平均时间"+num2str(timeSum/i)])

有约束的非线性规划问题

n=100000;
x_2=unifrnd(10,20,n,1);
x_1=x_2+10;
x_3=unifrnd(-10,16,n,1);     %通过两个不等式推出来的范围
res=-inf;
resX=0;
for i = 1:n
    X=[x_1(i),x_2(i),x_3(i)];
    if (-X(1)+2*X(2)+2*X(3)>=0) & (X(1)+2*X(2)+2*X(3)<=72)
        if X(1)*X(2)*X(3)>res
            res=X(1)*X(2)*X(3);
            resX=X;
        end
    end
end
disp(res);
disp(resX)

书店买书

min_money = +Inf;  % 初始化最小的花费为无穷大，后续只要找到比它小的就更新
min_result = randi([1, 6],1,5);  % 初始化五本书都在哪一家书店购买，后续我们不断对其更新
%若min_result = [5 3 6 2 3]，则解释为：第1本书在第5家店买，第2本书在第3家店买，第3本书在第6家店买，第4本书在第2家店买，第5本书在第3家店买  
n = 100000;  % 蒙特卡罗模拟的次数
M = [18	 39	29	48	59
        24	45	23	54	44
        22	45	23	53	53
        28	47	17	57	47
        24	42	24	47	59
        27	48	20	55	53];  % m_ij  第j本书在第i家店的售价
freight = [10 15 15 10 10 15];  % 第i家店的运费
for k = 1:n  % 开始循环
    result = randi([1, 6],1,5); % 在1-6这些整数中随机抽取一个1*5的向量，表示这五本书分别在哪家书店购买
    index = unique(result);  % 在哪些商店购买了商品，因为我们等下要计算运费
    money = sum(freight(index)); % 计算买书花费的运费
    % 计算总花费：刚刚计算出来的运费 + 五本书的售价
    for i = 1:5   
        money = money + M(result(i),i);  
    end
    if money < min_money  % 判断刚刚随机生成的这组数据的花费是否小于最小花费，如果小于的话
        min_money = money  % 我们更新最小的花费
        min_result = result % 用这组数据更新最小花费的结果
    end
end
min_money   % 18+39+48+17+47+20
min_result

导弹追踪离散化

v=200; % 任意给定B船的速度（后期我们可以再改的）
dt=0.0000001; % 定义时间间隔
x=[0,20]; % 定义导弹和B船的横坐标分别为x(1)和x(2)
y=[0,0]; % 定义导弹和B船的纵坐标分别为y(1)和y(2)
t=0; % 初始化导弹击落B船的时间
d=0; % 初始化导弹飞行的距离
m=sqrt(2)/2;   % 将sqrt(2)/2定义为一个常量，使后面看起来很简洁
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 导弹与B船的距离
for i=1:2
    plot(x(i),y(i),'.k','MarkerSize',1);  % 画出导弹和B船所在的坐标，点的大小为1，颜色为黑色(k)，用小点表示
    grid on;  % 打开网格线
    hold on;  % 不关闭图形，继续画图
end
axis([0 30 0 10])  % 固定x轴的范围为0-30  固定y轴的范围为0-10
k = 0;  % 引入一个变量  为了控制画图的速度（因为Matlab中画图的速度超级慢）
while(dd>=0.001)  % 只要两者的距离足够大，就一直循环下去。（两者距离足够小时表示导弹击中，这里的临界值要结合dt来取，否则可能导致错过交界处的情况）
    t=t+dt; % 更新导弹击落B船的时间
    d=d+3*v*dt; % 更新导弹飞行的距离
    x(2)=20+t*v*m;  y(2)=t*v*m;   % 计算新的B船的位置 （注：m=sqrt(2)/2）
    dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);  % 更新导弹与B船的距离
    tan_alpha=(y(2)-y(1))/(x(2)-x(1));   % 计算斜率，即tan(α)
    cos_alpha=sqrt(1/(1+tan_alpha^2));   % 利用公式：sec(α)^2 = (1+tan(α)^2)  计算出cos(α)
    sin_alpha=sqrt(1-cos_alpha^2);  % 利用公式： sin(α)^2 +cos(α)^2 = 1  计算出sin(α)
    x(1)=x(1)+3*v*dt*cos_alpha;   y(1)=y(1)+3*v*dt*sin_alpha;   % 计算新的导弹的位置
    k = k +1 ;  
    if mod(k,500) == 0   % 每刷新500次时间就画出下一个导弹和B船所在的坐标  mod(m,n)表示求m/n的余数
        for i=1:2
            plot(x(i),y(i),'.k','MarkerSize',1);
            hold on; % 不关闭图形，继续画图
        end
        pause(0.001);  % 暂停0.001s后再继续下面的操作
    end
    if d>50  % 导弹的有效射程为50个单位
        disp('导弹没有击中B船');
        break;  % 退出循环
    end
    if d<=50 & dd<0.001   % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）
        disp(['导弹飞行',num2str(d),'个单位后击中B船'])
        disp(['导弹飞行的时间为',num2str(t*60),'分钟'])
    end
end