MATLAB中的LogLog图(估计PI，误差与N)

Question

任务:使用Monte_Carlo方法，计算N=100,200,500,1000,2000,5000,10000,100000的PI近似值，并在对数图上绘制近似值与N的误差。其中Error=abs(实际值-近似值)。另外，使用另外两种无穷级数方法计算PI，并为N=10,20,50,100,200,500,1000,2000,5000,10000计算pi。在同一张图上绘制出所有2个公式和蒙特卡罗方法的对数图上的近似误差。

Estimating PI using M_C Method.
clear all
clc
close all
%n = linspace(0, 100000, 100)
n = [100, 200, 500, 1000, 2000, 5000, 10000, 100000]


c = 0;
t = 0;

for q=1:length(n)
    x = rand([1 n(q)]);
    y = rand([1 n(q)]);

    for i= 1:n(q)
        t = t+1;
        if x(i)^2 + y(i)^2 <= 1
            c = c+1;
            figure(2)
            k(i) = x(i);
            r(i) = y(i);
            hold on;
        else 
            figure(2)
            p(i) = x(i);
            j(i) = y(i);

        end
end
figure(1)
hold on
if n(q) == 1000

    plot(k, r, 'b.');
    plot(p, j, 'r.');
end
ratio = c/t;
PI= 4*ratio
error = abs(pi - PI)
figure(2)
loglog(n(q), error, '-b');




end

loglog(n, error, 's-')
grid on

%% Calculating PI using the James Gregory

%n = 10000;
n = [10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000]
d = 1;
c = 4;
for j = 1:n
    d = d + 2;
    c = c + ((-1)^(j))*(4)*(1/d);
end

PI_1 = c;
error = abs(n - PI_1);
loglog(n,error, '-s')
display(c);



%% Calculating PI using the Leibniz Formula for PI

%n = 10000;
n = [10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000]
d = 1;
c = 1;

for k = 1:n
    d = d + 2;
    c = c + ((-1)^k)*(1/d);
end

PI_2 = c*4;
error = abs(n - PI_2);
figure(3)
loglog(n, error, '-s')

我遇到的问题是loglog图没有显示预期的信息。

Answer 1

绘图讨论

对于蒙特卡罗绘图，这条线

loglog(n, error, 's-')

在for之后，-loop将覆盖由完成的所有绘图

loglog(n(q), error, '-b');

因为从来没有为figure(2)颁发过hold('on')。此外，在这两种情况下，由于样式选项和error不是向量的事实，绘图看起来会很奇怪：

• The call loglog(n, error, 's-')将生成一系列断开连接的框，因为n是一个向量，而error是一个标量；n将n的元素解释为独立的数据集，每个数据集都与相同的标量值error相关联(来自error for-loop的最后一次迭代的代码；<代码>d15是另一个称为<代码>d18的代码>具有类似的问题。由于样式需要一条“蓝色实线”，但每次迭代传递给loglog的数据是一个标量-标量对，因此不会出现任何内容。Matlab不能形成标量输入的直线(考虑线图plot(1,1,'-b')和圆图plot(1,1,'ob')作为另一个例子)。

这些问题可以通过将error更改为length(n)的向量来解决

error = zeros(1,length(n)); % before for-loop
...
    error(q) = abs(pi - PI); % inside the q-for-loop

并且仅在for-loop之后执行loglog绘图(这也是一个性能提升，因为绘图调用相对于计算而言是繁重的)。

性能讨论

说到性能(为了加快蒙特卡罗的速度)，蒙特卡罗集成的最大优点，除了不会屈从于curse of dimensionality之外，还在于它的可笑的并行化(即，可矢量化)性质。

这是一件很好的事情，因为香草蒙特卡洛需要大量的样本才能得到准确的结果。此外，Matlab的logical indexing允许使用一种很好的语义方法来提取满足多个条件的值。

也就是说，可以使用以下代码对蒙特卡罗方法的i-for-loop进行审查：

% % ----- i-for-loop replacement

%   Determine location of points
inCircle = (x.^2 + y.^2) <= 1;

%   k = xIn, r = yIn
xIn = x(inCircle);
yIn = y(inCircle);

%   p = xOut, j = yOut
xOut = x(~inCircle);    % or x(not(inCircle));
yOut = y(~inCircle);    % or y(not(inCircle));

% % ----- end of i-for-loop replacement

%   Calculate MC pi and error
ratio    = nnz(inCircle)/n(q);
PI       = 4*ratio;
error(q) = abs(pi - PI);