% Biased vs unbiased variance estimators
% Textbook: Part I, "Descriptive Statistics", exercise "simulation: both variance estimators"
%
% k experiments, each of n samples from N(0,1). In every experiment compute
%   s2_biased   = sum((x - mean(x)).^2)/n       -> E = (n-1)/n * sigma^2
%   s2_unbiased = sum((x - mean(x)).^2)/(n-1)   -> E = sigma^2
% MSE about sigma^2: unbiased 2/(n-1)*sigma^4, biased (2n-1)/n^2*sigma^4.
% The unbiased estimator is wider, so its MSE is the larger of the two.

clearvars; close all
rng(0)

k = 1e4;        % number of experiments
n = 10;         % samples per experiment
sigma2 = 1;     % true variance

x = sqrt(sigma2)*randn(n, k);       % one experiment per column
s2_b = var(x, 1);                   % divide by n
s2_u = var(x);                      % divide by n-1 (MATLAB default)

fprintf('%-10s %12s %12s %12s %12s\n', '', 'mean, sim', 'mean, theory', 'MSE, sim', 'MSE, theory')
fprintf('%-10s %12.4f %12.4f %12.4f %12.4f\n', 'biased', ...
    mean(s2_b), (n-1)/n*sigma2, mean((s2_b - sigma2).^2), (2*n-1)/n^2*sigma2^2)
fprintf('%-10s %12.4f %12.4f %12.4f %12.4f\n', 'unbiased', ...
    mean(s2_u), sigma2, mean((s2_u - sigma2).^2), 2/(n-1)*sigma2^2)

figure
histogram(s2_b, 'Normalization', 'pdf', 'BinWidth', 0.1)
hold on
histogram(s2_u, 'Normalization', 'pdf', 'BinWidth', 0.1)
xline(sigma2, 'k--', 'LineWidth', 1.5)
hold off
xlabel('variance estimate')
ylabel('pdf')
legend('biased, /n', 'unbiased, /(n-1)', 'true \sigma^2')
grid on