Экспериментальное вычисление "пи" методом Монте-Карло

a = 0;
b = 1;
n = 1e6;
x = a + (b-a)*rand(1,n);
y = a + (b-a)*rand(1,n);
fun = sqrt(1-x.^2);

ns = 0; % счетчик попаданий внутрь окружности
for i=1:n
    if y(i) <= fun(i)
        ns = ns +1;
    end
end

disp(4*ns/n)

n = 1e6;
x = rand(1,n);
fun = sqrt(1-x.^2);
y = rand(1,n);
ns = sum(y <= fun); % счетчик попаданий внутрь окружности

plot(x,fun,'.',x(y <= fun),y(y <= fun),'g.',x(y > fun),y(y > fun),'r.'),axis('square')