% Moments of Y = X^2 + 1 for an exponential X, as a function of N
% Textbook: Part I, "Continuous Random Variable", exercise (2022, Classwork)
%
% X ~ Exp(lambda), lambda = 1/2. Analytically,
%   E[Y] = 2/lambda^2 + 1 = 9,   Var[Y] = 20/lambda^4 = 320.
% The simulated values approach the theory as N grows, but slowly: Y has
% heavy tails, so even N = 1e6 is not exact.

clearvars; close all
rng(0)

lambda = 1/2;
N = 1e6;
x = -log(rand(1, N))/lambda;    % Exp(lambda) by inverse-CDF (or exprnd(1/lambda,1,N))
y = x.^2 + 1;

M = round(logspace(2, 6, 20));  % growing number of experiments
E = zeros(size(M));
V = zeros(size(M));
for k = 1:length(M)
    E(k) = mean(y(1:M(k)));
    V(k) = var(y(1:M(k)));
end

figure
tiledlayout(1, 3)
nexttile
semilogx(M, E, 'o-', M, 9*ones(size(M)), 'r--', 'LineWidth', 1.2)
xlabel('N'); title('E[Y]'); legend('simulation', 'theory'); grid on
nexttile
semilogx(M, V, 'o-', M, 320*ones(size(M)), 'r--', 'LineWidth', 1.2)
xlabel('N'); title('Var[Y]'); legend('simulation', 'theory'); grid on
nexttile
histogram(y, 'Normalization', 'pdf', 'BinLimits', [1 50])
xlabel('y'); title('f_Y(y)'); grid on

fprintf('E[Y]:   theory = 9,   simulation (N = %g) = %.3f\n', N, E(end))
fprintf('Var[Y]: theory = 320, simulation (N = %g) = %.3f\n', N, V(end))