% Fitting a distribution to a pdf histogram
% Textbook: Part I, "Statistical Plots", pdf histograms
%
% Estimate the parameters of a model distribution from the data, then
% overlay the model density on the pdf histogram:
%   exponential: lambda_hat = 1/mean(x)
%   Gaussian:    mu_hat = mean(x), sigma_hat = std(x)
% (fitdist in the Statistics Toolbox does the same.)

clearvars; close all
rng(0)

N = 1e4;

% Exponential data, lambda = 0.5
xe = -log(rand(1, N))/0.5;
lambda_hat = 1/mean(xe);

% Gaussian data, mu = 3, sigma = 2
xg = 3 + 2*randn(1, N);
mu_hat = mean(xg);
sigma_hat = std(xg);

fprintf('exponential: lambda = 0.5, estimate = %.4f\n', lambda_hat)
fprintf('Gaussian:    mu = 3, estimate = %.4f; sigma = 2, estimate = %.4f\n', mu_hat, sigma_hat)

figure
tiledlayout(1, 2)
nexttile
histogram(xe, 'Normalization', 'pdf')
hold on
t = linspace(0, max(xe), 200);
plot(t, lambda_hat*exp(-lambda_hat*t), 'r', 'LineWidth', 1.5)
hold off
xlim([0 15])
xlabel('x'); ylabel('f_X(x)')
title(sprintf('Exp, \\lambda = %.3f', lambda_hat))
grid on

nexttile
histogram(xg, 'Normalization', 'pdf')
hold on
t = linspace(min(xg), max(xg), 200);
plot(t, exp(-(t - mu_hat).^2/(2*sigma_hat^2))/(sqrt(2*pi)*sigma_hat), 'r', 'LineWidth', 1.5)
hold off
xlabel('x'); ylabel('f_X(x)')
title(sprintf('N(%.2f, %.2f^2)', mu_hat, sigma_hat))
grid on