main =>
   _ = random2(), % random seed
   G = [gaussian_dist(1,0.5) : _ in 1..1000],
   println(first_10=G[1..10]),
   println([mean=avg(G),stdev=stdev(G)]),
   nl.

% Gaussian (Normal) distribution, Box-Muller algorithm
gaussian01() = Y =>
    U = frand(0,1),
    V = frand(0,1),
    Y = sqrt(-2*log(U))*sin(2*math.pi*V).

gaussian_dist(Mean,Stdev) = Mean + (gaussian01() * Stdev).

% Variance of Xs
variance(Xs) = Variance =>
    Mu = avg(Xs),
    N  = Xs.len,
    Variance = sum([ (X-Mu)**2 : X in Xs ]) / N.

% Standard deviation
stdev(Xs) = sqrt(variance(Xs)).