%% （2）代码部分 l = 0.520; % 针的长度（任意给的） a = 1.314; % 平行线的宽度(大于针的长度l即可) n = 1000000; % 做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， 并打开图形的边框 fori=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)])
%% （3） 由于一次模拟的结果具有偶然性，因此我们可以重复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; fori=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)])
1
rand(5,4)
这行代码生成一个大小为 5x4 的矩阵，其中的元素是在区间 [0, 1]
内均匀分布的随机数。
1 2
matlabCopy code rand(3)
这行代码生成一个大小为 3x3 的方阵，其中的元素同样是在区间 [0, 1]
内均匀分布的随机数。
1 2
matlabCopy code -2 + rand(3,2) * (2 - (-2))
这行代码生成一个大小为 3x2 的矩阵，其中的元素是在区间 [-2, 2]
内均匀分布的随机数。
1 2
matlabCopy code unifrnd(-2,2,3,2)
这是另一种生成大小为 3x2 的矩阵，元素在区间 [-2, 2]
内均匀分布的方法。
现在，我们来看实际的代码部分：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
matlabCopy codel = 0.520; % 针的长度 a = 1.314; % 平行线的宽度（应大于针的长度） n = 1000000; % 投针实验的次数 m = 0; % 记录针与平行线相交的次数
x = rand(1, n) * a / 2; % 在 [0, a/2] 内生成 n 个随机数，表示针的中点距离最近平行线的距离 phi = rand(1, n) * pi; % 在 [0, pi] 内生成 n 个随机数，表示针与最近平行线的夹角
% 循环，检查每次实验中针是否与平行线相交 fori = 1:n if x(i) <= l / 2 * sin(phi(i)) m = m + 1; % 如果相交，计数加一 end end
p = m / n; % 针与平行线相交的概率 mypi = (2 * l) / (a * p); % 使用蒙特卡罗方法估算 pi
disp(['蒙特卡罗方法得到的 pi 估计值：', num2str(mypi)])
该代码模拟了布丰投针实验，其中一根长度为 l
的针被投在两条间距为 a
的平行线上。通过估算针与线相交的概率，利用蒙特卡罗方法估计 π 的值。
%% （2）代码部分（在成功的条件下的概率） n = 100000; % n代表蒙特卡罗模拟重复次数 a = 0; % a表示不改变主意时能赢得汽车的次数 b = 0; % b表示改变主意时能赢得汽车的次数 fori= 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表示没有获奖的次数 fori= 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)]);
figure(1) % 新建一个编号为1的图形窗口 plot(coord(:,1),coord(:,2),'o'); % 画出城市的分布散点图 fori = 1:n text(coord(i,1)+0.01,coord(i,2)+0.01,num2str(i)) % 在图上标上城市的编号（加上0.01表示把文字的标记往右上方偏移一点） end hold on % 等一下要接着在这个图形上画图的
