MATLAB実装!! 確率使って正方形と扇から円周率を作る

2021/6/18 コード更新 (MATLAB R2021aで実装)

「大学生になってパソコンを持ったら何でもできるようになる」って思ってたんだけど, なんだかんだ"離散"という制約により, できないことも多い気がする...

円周率と確率

円周率をパソコンでもとっめる方法はたくさんあるが, 正方形と円に一様な点を打っていく方法はとても単純で有名である. 今回は, 以下のような正方形と1/4円を使って円周率を求める.

Code

% Create circular sector
t = 0:0.01:pi/2;
xc = cos(t);
yc = sin(t);

% Plot circular sector and rectangle
figure
hold on
plot([0,1,1,0,0],[0,0,1,1,0],'b-','LineWidth',1.5);
plot(xc,yc,'k-','LineWidth',1.5);
hold off
grid on
axis equal
title('Field');
xlabel('x');
ylabel('y');
xlim([-0.25 1.25]);
ylim([-0.25 1.25]);

上の図において, 青色の正方形の面積は,
$$ A_c = 1\times 1 = 1 $$
であり, 黒色の扇形の面積は
$$ A_f = \frac{1}{4} \times1\times 1\times \pi = \frac{\pi}{4} $$
となる. ここに, $[0,1]$の範囲で一様分布な点を打っていくと, それが扇形の範囲の中に入る確率$r$は,
$$ r = \frac{A_f}{A_c}=\frac{\pi}{4} $$
となる. つまり, この範囲に点を打っていき, それが扇型の中に入っている確率の4倍が円周率となる.

MATLABで実装

これをMATLABで実装する.

乱数とデータの作成

MATLABにはたくさんの乱数発生器が用意されているが, 今回は, rand関数$0 \leq x,y \leq 1$の点を作成する.

Code

%% Setting Part
mi = 500;      % Max Iteration Number
error = 0.001; % Acceptable Error
stop = 2;      % Howmany times error should be under mi

p = [];    % data
s = 0;     % Distance from Origin (0,0)
r = 0;     % ratio*4
judge = 0; % Count allowable error

% Iteration Part
i = 1;
while (judge < stop && i < mi)
    p = [p,rand(2,1)];           % Generate data
    s(i) = sqrt(sum(p(:,i).^2)); % Caluculate Distance from Origin
    inner = s(s <= 1);           % Judge inside or outside of circular sector
    r(i) = (numel(inner)/i)*4;   % Ratio*4 shoud be pi
    
    % Calculate error
    if (abs(r(i)-pi) < error) 
        judge = judge + 1; 
    end
    i = i + 1;
end
count = i-1;

可視化

上のコードの結果をグラフにします.

Code

%% Visualize Part
% Scatter Plot
figure;
subplot(5,1,[1,2,3]);
hold on
plot([0,1,1,0,0],[0,0,1,1,0],'b-','LineWidth',1.5);
plot(xc,yc,'k-','LineWidth',1.5);
p1 = plot(p(1,1),p(2,1),'g.','MarkerSize',8);
hold off
axis equal
grid on
title('Scatter Plots');
xlabel('x');
ylabel('y');

% Ratio*4 plot
subplot(5,1,4);
hold on
p2 = plot(r(1),'g-','LineWidth',1.1);
yline(pi,'r-','LineWidth',0.7);
hold off
grid on
xlim([1 count]);
ylim([0 5]);
xlabel('Iteration Number');
ylabel('Ratio (\pi)');
title('Result');

% Error from pi plot
subplot(5,1,5);
p3 = semilogy(abs(r(1)-pi),'r-','LineWidth',1.1);
grid on
xlim([1 count]);
ylim([min(abs(r-pi)) 3]);
xlabel('Iteration Number');
ylabel('Error');
title('Error');

% Create Animation
for i=1:count
    p1.XData = p(1,1:i);
    p1.YData = p(2,1:i);
    p2.YData = r(1:i);
    p3.YData = abs(r(1:i)-pi);
    drawnow
end

結果, 以下のように円周率に向かって誤差が小さくなっていきます.

円周率を正方形と扇形から求めよう!!

ここで使った乱数発生器はMATLAB rand関数. 他にもMATLABにはいろんな種類の乱数発生器があって, とっても便利!!

ちなみに私一押しはrandn関数.
測定誤差を含んだデータを簡単に再現できます (後日のネタにとっておきます).#MATLAB #uec #電通大 #Qiita pic.twitter.com/oz1vvTNqfB

— ヒビヤ@MATLAB Ambassador (@hibs_MATLAB_Amb) June 17, 2021