蒙特卡罗简单学习

Posted 2021-09-15 csp_

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了蒙特卡罗简单学习相关的知识，希望对你有一定的参考价值。

布丰投针实验求π

代码：

l =  0.520;   % 针长
a = 1.314;    % 平行线的宽度(大于针长即可)
n = 10000;    % 做n次投针试验
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;  
        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)])

结果：

蒙特卡罗方法得到pi为：3.1358

可以重复100次后再来求一个平均的pi：

result = zeros(100,1);  % 初始化保存100次结果的矩阵
l =  0.520;     a = 1.314;
n = 1000000;    
for num = 1:100
    m = 0;  
    x = rand(1, n) * a / 2 ;
    phi = rand(1, n) * pi;
    for i=1:n
        if x(i) <= l / 2 * sin(phi (i))
            m = m + 1;
        end
    end
    p = m / n;
    mypi = (2 * l) / (a * p);
    result(num) = mypi;  % 把求出来的myphi保存到结果矩阵中
end
mymeanpi = mean(result);  % 计算result矩阵中保存的100次结果的均值
disp(['蒙特卡罗方法得到pi为：', num2str(mymeanpi)])

结果：

蒙特卡罗方法得到pi为：3.1412

错排问题求e

知识：

randperm(5)  % 生成1-5组成的一个随机序列
%      3     5     1     2     4
%      1     4     5     3     2

% 假设a是一个向量，那么find(a)可以用来返回这个向量中非零元素的下标，如果a中所有元素都为0，则返回空值
find([1,5,6,0,8,0,-5])  %      1     2     3     5     7
find([0,0,0,0,0])  %   空的 1×0 double 行矢量

% (3) 矩阵（或向量）和常量的比较运算可返回逻辑矩阵（或向量）(元素全为0和1)
[1,5,6,0,8,0,-5] > 0      %    1   1   1   0   1   0   0
[1,5,6,0,8,0,-5] == 0    %    0   0   0   1   0   1   0

% (4) isempty(A)函数可以用来判断A是否为空, 如果A为空， isempty(A) 返回逻辑值1(true)，否则返回逻辑值0(false)。
isempty(find([0,0,0,0,0]))   %    1
isempty(find([0,1,0,0,0]))   %    0
isempty([0,0,0,0,0])  %它不是空矩阵（空矩阵是指里面没有元素）

求解：

clear;clc
tic  %计算tic和toc中间部分的代码的运行时间
n = 1000000; 
m = 0;   % 每个人拿到的都不是自己卡片的次数（频数）
people = 100;   % 假设一共有100个人玩这个游戏 
for i = 1: n  % 开始循环
    if isempty(find(randperm(people) - [1:people] == 0))  % 如果每个人拿到的都不是自己的卡片
        m = m + 1;  % 那么次数就加1
    end
end
frequency = m / n;  % 每个人拿到的都不是自己卡片的频率（概率）
disp(['自然常数e的蒙特卡罗模拟值为：', num2str(1 / frequency)])  % 注：自然常数真实值约为2.7182
toc

结果：

自然常数e的蒙特卡罗模拟值为：2.7167
历时 5.235055 秒

三门问题

蒙特卡罗求在成功的条件下的概率

n = 100000;  % 蒙特卡罗模拟重复次数
a = 0;  % 不改变主意时能赢得汽车的次数
b = 0;  % 改变主意时能赢得汽车的次数
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)]);

输出：

蒙特卡罗方法得到的不改变主意时的获奖概率为：0.33452
蒙特卡罗方法得到的改变主意时的获奖概率为：0.66548

蒙特卡罗考虑失败情况的概率(无条件概率)

n = 100000;  % 蒙特卡罗模拟重复次数
a = 0;  % 不改变主意时能赢得汽车的次数
b = 0;  % 改变主意时能赢得汽车的次数
c = 0;  % 没有获奖的次数
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)]);

结果：

蒙特卡罗方法得到的不改变主意时的获奖概率为：0.16723
蒙特卡罗方法得到的改变主意时的获奖概率为：0.33271
蒙特卡罗方法得到的没有获奖的概率为：0.50006

模拟排队问题

问题1的代码：

tic  %计算tic和toc中间部分的代码的运行时间
i = 1;  
w = 0;  % w用来表示所有客户等待的总时间，初始化为0
e0 = 0;  c0 = 0;   
x(1) = exprnd(10);  
c(1) = c0 + x(1);  
b(1) = c(1); 
while b(i) <= 480  % 开始设置循环，只要第i个顾客开始服务的时间(时刻)小于480，就可以对其服务（银行每天工作8小时，折换为分钟就是480分钟）
    y(i) = normrnd(10,2); % 第i个客户的服务持续时间，服从均值为10方差为4(标准差为2)的正态分布
    if y(i) < 1  
        y(i) = 1;
    end
    e(i) = b(i) + y(i); % 第i个客户结束服务的时间 = 第i个客户开始服务的时间 + 第i个客户的服务持续时间
    wait(i) = b(i) - c(i); % 第i个客户等待的时间 = 第i个客户开始服务的时间 - 第i个客户到达银行的时间
    w = w + wait(i); % 更新所有客户等待的总时间
    i = i + 1; % 增加一名新的客户
    x(i) = exprnd(10); % 这位新客户和上一个客户到达的时间间隔
    c(i) = c(i-1) + x(i); % 这位新客户到达银行的时间 = 上一个客户到达银行的时间 + 这位新客户和上一个客户到达的时间间隔
    b(i) = max(c(i),e(i-1)); % 这个新客户开始服务的时间取决于其到达时间和上一个客户结束服务的时间
end
n = i-1; 
t = w/n; % 客户的平均等待时间
disp(['银行一天8小时一共服务的客户人数为: ',num2str(n)])
disp(['客户的平均等待时间为: ',num2str(t)])
toc  %计算代码运行时间

问题2的代码：

tic  
day = 100;  % 假设模拟100天
n = zeros(day,1); % 每日接待客户数
t = zeros(day,1); % 每日客户平均等待时长
for k = 1:day
    i = 1;  
    w = 0;  
    e0 = 0;  c0 = 0;   
    x(1) = exprnd(10);  
    c(1) = c0 + x(1);  
    b(1) = c(1); 
    while b(i) <= 480  
        y(i) = normrnd(10,2); 
        if y(i) < 1  
            y(i) = 1;
        end
        e(i) = b(i) + y(i); 
        wait(i) = b(i) - c(i); 
        w = w + wait(i); 
        i = i + 1; % 
        x(i) = exprnd(10); 
        c(i) = c(i-1) + x(i); 
        b(i) = max(c(i),e(i-1)); 
    end
    n(k) = i-1; % n(k)表示银行第k天服务的客户人数
    t(k) = w/n(k); % t(k)表示该银行第k天客户的平均等待时间
end
disp([num2str(day),'个工作日中，银行每日平均服务的客户人数为: ',num2str(mean(n))])
disp([num2str(day),'个工作日中，银行每日客户的平均等待时间为: ',num2str(mean(t))])
toc  %计算代码运行时间

解有约束的非线性规划问题

求满足以下约束条件下函数 $f(x)=x_{1}x_{2}x_{3}$ 的最大值
$t.\\begin{cases}-x_{1}+2x_{2}+2x_{3}>=0\\\\ x_{1}+x_{2}+2x_{3}<=72\\\\ 10<=x_{2}<=20\\\\ x_{1}-x_{2}=10\\end{cases}$

五子棋对弈——MCTS学习

强化学习读书手札：动态规划(DP)&蒙特卡洛(MC)&时序差分(TD)区别

IOS开发－OC学习－常用功能代码片段整理

theano学习指南--混合蒙特卡洛采样(翻译)

蒙特卡罗方法python简单实现

与 openmp 的 10 维蒙特卡罗集成

蒙特卡罗简单学习

目录

布丰投针实验求π

错排问题求e

三门问题

模拟排队问题

解有约束的非线性规划问题