clc
clear 
close all

N=1:10:500; %variare il numero di dati
r = 0.5;
iterations = 100; %variare il numero di iterazioni

pigreco_tmp = ones(iterations,1);
pigreco_stima = ones(length(N),1);
percentili = ones(length(N),2);

for i = 1:length(N)
    for j=1:iterations
        
        %estrazione punti (x,y) scelti in modo equiprobabile nel quadrato 
        px=rand(N(i),1);
        py=rand(N(i),1);
        
        f=@(x,y) (x-r).^2 + (y-r).^2;
        F=f(px,py);
        
        % punti nel cerchio
        N_cerchio=length(find(F<=r^2));
        
        % A_cerchio/A_quadrato = N_cerchio/N_tot
        pigreco_tmp(j)=N_cerchio/(r^2*N(i));
    end
    
    pigreco_stima(i) = median(pigreco_tmp);
    percentili(i,:)= prctile(pigreco_tmp,[2.5 97.5]);
end

plot(N,pigreco_stima,'k')
hold 
plot(N,pi*ones(size(N)),'r')
hold on
plot(N,percentili(:,1),'--k')
hold on
plot(N,percentili(:,2),'--k')

grid
xlabel('N')
ylabel('pi')
legend('pi estimation','pi true','95% limits')