% Monte-Carlo error vs the number of experiments k
% Textbook: Part I, "Descriptive Statistics", Monte-Carlo simulation
%
% A statistic g with E[g] = 0 and std sigma = 1 is simulated k times.
% - The simulated mean fluctuates around E[g] within about sigma/sqrt(k).
% - The relative error of the simulated std decays as 1/sqrt(2k).
% Ten independent simulation runs are shown; each extra digit of accuracy
% costs 100 times more experiments.

clearvars; close all
rng(0)

n_runs = 10;
k_max = 1e6;
g = randn(n_runs, k_max);           % g ~ N(0,1), one run per row

k = 1:k_max;
m = cumsum(g, 2)./k;                % running mean
s = sqrt(cumsum(g.^2, 2)./k - m.^2);% running std

idx = unique(round(logspace(1, log10(k_max), 80)));

figure
tiledlayout(1, 2)
nexttile
semilogx(idx, m(:, idx), 'Color', [0.5 0.5 0.8])
hold on
semilogx(idx, 1./sqrt(idx), 'k--', idx, -1./sqrt(idx), 'k--', 'LineWidth', 1.5)
hold off
xlabel('k')
title('simulated mean, envelope \pm1/\surd k')
grid on

nexttile
loglog(idx, abs(s(:, idx) - 1), 'Color', [0.5 0.5 0.8])
hold on
loglog(idx, 1./sqrt(2*idx), 'k--', 'LineWidth', 1.5)
hold off
xlabel('k')
title('relative std error, 1/\surd(2k)')
grid on

fprintf('k = %g: std of simulated mean over runs = %.2e (theory %.2e)\n', ...
    k_max, std(m(:, end)), 1/sqrt(k_max))