モンテカルロ法によるPIの推定(ビフォンの問題)

%ビフォンの問題
No = 100;
D = 0.2;
L = 0.1;
line = L/2:D:1-L/2
for j=1:No
N = j*100;
y = rand(N,1);
theta = rand(N,1)*pi/2;
y1 = y + L/2*cos(theta);
y2 = y - L/2*cos(theta);
s = 0;
for i=1:length(line)
s = s + sum( (line(i) >= y1) == (line(i) <= y2) );
end
%
% p = 2*L/(pi*D)
%
p = s/N
PI(j) = 2*L/(D*p);
NS(j) = N;
end

hold on;
plot(NS,PI,'.k');
plot(NS,ones(No,1)*pi,'--r');
title('ビフォンの問題');
xlabel('繰り返し回数');
ylabel('\piの推定');
grid on;
hold off;
print('-dpng','-r80','pi_buffo.png');