function MonteCarloExample1 (A, N, M)
binsize = 0.015*A;
b = 0.1*A:binsize:3*A;
close all;
for k = 1:length(N)
    Aest = zeros(1,M);
    for m = 1:M
        w = sqrt(A)*randn(1,N(k));
        x = A + w;
        Aest(m) = -0.5+sqrt(mean(x.^2)+0.25);
    end;
    meanAest(k) = mean(Aest);
    varAest(k) = var(Aest,1);
    varAtheoretical(k) = (A.^2)/(N(k)*(A+0.5));
    subplot(3,ceil(length(N)/3),k);
    s = sprintf('N = %d, M = %d',N(k),M);
    title(s);hold on;
    h = hist(Aest,b);
    plot(b,(h/M)/binsize,'b-*');
    p=theoreticalPDF(A,N(k),b);
    plot(b,p,'r-*');
end;
figure;bar(meanAest);title('Mean of estimator');
figure;bar([varAest; varAtheoretical]');title('Variance of estimator');legend('Empirical','Theoretical');
figure;bar([varAest.*N; varAtheoretical.*N]');title('N times variance of estimator');legend('Empirical','Theoretical');
end

function p = theoreticalPDF(A,N,x)
var = (A.^2)/(N*(A+0.5));
p = (1/sqrt(2*pi*var))*exp((-0.5/var)*((x-A).^2));
end