new file mode 100644
--- /dev/null
+++ b/drawmedcouple.m
@@ -0,0 +1,38 @@
+
+# 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