基于MATLAB，使用蒙特卡洛法进行Pi值的估算

Pi值的估算

估算方法：蒙特卡洛法

实验原理：

估算模型：

在边长为2的正方形内随机生成点，用随机掉落在半径为1的圆内点数，除以掉落在其外接正方形内点数，来近似估算pi。

以下是随机生成点数为10000,20000,30000。其pi的估值分别为3.1872 ，3.1354 ，3.1455。

实验次数N与绝对误差的关系（近似为正态分布）

一次随机实验中，点落在圆内的概率为p=pi/4。n为实验次数，置信度至少为95%（pi_alpha=0.05），则绝对误差公式：delta=sqrt((2*p*(1-p))/(n*erf(1-alpha)^2))以，实验次数从10000到10000000每次实验次数增加10000次时，绝对误差和实验次数的图像。

这个图像表明随实验次数N求得的pi的均值，符合正态分布。

以下该图像为实验次数N增大，同时设定样本数也为N，发现其也符合正态分布，猜测原因是，在实验次数增加时，p趋近真实值，符合大数定理，则p接近Ep，与求样本的均值分布都满足正态分布。

完整代码

clear all
close all
%参数
N=10000:10000:10000000;
numbber=1000;
numbberx=1:numbber;
N_length=length(N);
Circle_X=1;
Circle_Y=1;
Circle_R=1;
Rectangle_X=2;
Rectangle_Y=2;
Estimate_pi_sample=[];
alpha=0.05;
p=pi/4;
delta_y=[];
Estimate_pi_sample_sum=[];
%画圆
for i=1:numbber
    theta=0:pi/20:2*pi;
    Trajectory_Circle_X=1+cos(theta);
    Trajectory_Circle_Y=1+sin(theta);
%     subplot(3,1,i);
%     plot(Trajectory_Circle_X,Trajectory_Circle_Y);
%     hold on;
    %绘制点
    Point_Incircle=0;
    for n=1:numbber
        Random_Point_X=2*rand;
        Random_Point_Y=2*rand;
        if (Random_Point_X-1)^2+(Random_Point_Y-1)^2<1
             Point_Incircle=Point_Incircle+1;
%             plot(Random_Point_X,Random_Point_Y,'.','Color',[0,0,1]);
%         else
%             plot(Random_Point_X,Random_Point_Y,'.','Color',[1,0,0]);
        end
    end
    %估算pi
    Estimate_pi=(Point_Incircle/numbber)*4;
    Estimate_pi_sample=[Estimate_pi_sample Estimate_pi];
    %pi均值
    Estimate_pi_sample_sum=[Estimate_pi_sample_sum sum(Estimate_pi_sample)/i];

end
%验证大数定理验证，当n增大时，当N越大时，Estimate_pi越接近pi，绝对误差近似正态分布
%验证大数定理验证，当n确定时，当N越大时，Estimate_pi越接近pi，绝对误差近似正态分布
for i=1:length(N)
    delta=sqrt((2*p*(1-p))/(N(i)*erfinv(1-alpha)^2));
    delta_y=[delta_y delta];
end
subplot(211)
plot(N,delta_y)
xlabel('实验次数N');
ylabel('理论绝对误差delta');
title('实验次数和理论绝对误差的关系');
subplot(212)
plot(numbberx,abs(Estimate_pi_sample_sum-pi))
xlabel('实验次数N');
ylabel('实际绝对误差delta');
title('实验次数和实际绝对误差的关系');