# Set the seed, to make this "random" sampling deterministic
randg("seed", 0);

z = gamrnd (3, 2, 5e3, 1);

# Plot the gamma values
subplot(2,1,1);
hist(z,20, "facecolor", [1, 0.5, 0.5]);
title("gamma distribution with alpha=3 and beta=2", "fontsize", 20);
set(gca, "fontsize", 14)

# Naive medcouple algorithm
n = length(z);
n2 = ceil(n/2);
z = sort(z, "descend");
zmed = median(z);
z -= zmed;
zplus = z(z >= 0);
zminus = z(z <= 0)';
zz = (zplus + zminus)./(zplus - zminus);
zz(isnan(zz)) = 0;
zzsort = sort(zz(:), "descend");
medc_idx = ceil(n2^2/2);
mc = zzsort(medc_idx);

# Plot the medcouple distribution
subplot(2,1,2);
hist(zzsort,20, "facecolor", [0.5, 0.5, 1]);

# Draw the medcouple line on the medcouple distribution
hold on;
plot([mc,mc], [0,5e5], "linewidth", 2.5, "color", [0.95, 0.95, 0]);
title("distribution of h(x_i, x_j) values", "fontsize", 20);
set(gca,"xtick", [-1, -0.5, mc, 0.5, 1], "fontsize", 14);
hold off

print medcouple-distribution.png