% Sample skewness and excess kurtosis
% Textbook: Part I, "Descriptive Statistics", skewness and tails
%
% skewness = E[(X-mu)^3]/sigma^3, excess kurtosis = E[(X-mu)^4]/sigma^4 - 3.
% Four distributions, the last three scaled to zero mean and unit variance:
%   log-normal (positive skew), Laplace (heavy tails, +3),
%   normal (0, 0), uniform (light tails, -6/5).

clearvars; close all
rng(0)

N = 1e6;
s_ln = 0.8;                                 % log-normal shape
data = { ...
    exp(1.2 + s_ln*randn(1, N)), ...        % log-normal
    (log(rand(1, N)) - log(rand(1, N)))/sqrt(2), ... % Laplace, b = 1/sqrt(2)
    randn(1, N), ...                        % normal
    sqrt(3)*(2*rand(1, N) - 1)};            % uniform on [-sqrt(3), sqrt(3)]
names = {'log-normal', 'Laplace', 'normal', 'uniform'};

w = exp(s_ln^2);
skew_th = [(w + 2)*sqrt(w - 1), 0, 0, 0];
kurt_th = [w^4 + 2*w^3 + 3*w^2 - 6, 3, 0, -6/5];

figure
tiledlayout(2, 2)
fprintf('%-11s %10s %10s %10s %10s\n', '', 'skew, sim', 'skew, th', 'kurt, sim', 'kurt, th')
for i = 1:4
    x = data{i};
    z = (x - mean(x))/std(x, 1);
    skew = mean(z.^3);              % skewness(x) in the Statistics Toolbox
    kurt = mean(z.^4) - 3;          % kurtosis(x) - 3
    fprintf('%-11s %10.3f %10.3f %10.3f %10.3f\n', names{i}, skew, skew_th(i), kurt, kurt_th(i))

    nexttile
    histogram(x, 200, 'Normalization', 'pdf', 'EdgeColor', 'none')
    hold on
    xline(mean(x), 'b', 'mean', 'LineWidth', 1.2)
    xline(median(x), 'r', 'median', 'LineWidth', 1.2)
    hold off
    title(sprintf('%s: skew %.2f, kurt %.2f', names{i}, skew, kurt))
    if i == 1, xlim([0 25]); else, xlim([-4 4]); end
    grid on
end