本文已参与「新人创作礼」活动，一起开启掘金创作之路。

一个奇奇怪怪的任务，关于Matlab 蒙特卡洛模拟

客户提出了奇奇怪怪的要求，于是我也写出了奇奇怪怪的代码，如果您对奇奇怪怪的要求没有兴趣，那么请忽略这篇奇奇怪怪的文章。

任务描述

这次的任务与其说是编程，其实更像是阅读理解，客户的要求如下：

• matlab计算

• 过程：生成1行100列的矩阵，每个元素只能是0或者1

  • 其中取1的概率为pn，0的概率为1-pn。
  • 第80列位置取1的概率是1-pn。重复50次这个过程，然后将50次结果累加。

• pn从0%，1%，2%，3%逐渐增加到100%，求出80这个位置数值最大的概率ps

• 画坐标图，横轴为pn,纵轴为ps。

实现代码

pn_ps = zeros(1,101);
time = 1000; % 实验次数，默认100
for pn=1:101
   for i = 1:time       %实验time次 
       out = zeros(1,100);
       for t = 1:50     % 重复50次这个过程
           for j=1:100  % 对于每1列 
               if rand>(pn-1)/100 % 取1的概率为pn rand生成随机数(0-1)
                   out(j)=out(j)+0;
               else
                   out(j)=out(j)+1; % 取1时累加
               end
           end
       end
       out(80) = 50-out(80); % 第80列结果取反
       if out(80)==max(out)  % 第80列如果数值最大，则统计
           pn_ps(pn) = pn_ps(pn)+1;
       end
   end
   pn_ps(pn)=pn_ps(pn)/time*100; % 概率归一化
end
plot((0:1:100),pn_ps/time*100)  %做图
xlabel('Pn') % x轴
ylabel('Ps') % y轴

结果

后记

这个结果有些类似于深度学习的激活函数程序上面的段落可以用随机向量函数rand(1,100)+判断向量的值大于特定值赋1，小于特定值赋0的方法替代，应该能降低时间复杂度。

for j=1:100  % 对于每1列 
 if rand>(pn-1)/100 % 取1的概率为pn rand生成随机数(0-1)
    out(j)=out(j)+0;
 else
    out(j)=out(j)+1; % 取1时累加
 end
end

time取1000时，运行时间已经有些长了，如果取10000，要运行很久才能获得较为平滑的图像，这是蒙特卡洛法的缺点，感觉客户不应该用这种方法，应该直接用公式计算概率。