/* 

  Random numbers (Rosetta code) in Picat.

  http://rosettacode.org/wiki/Random_numbers
  """
  Task

  Generate a collection filled with   1000   normally distributed random (or pseudo-random) 
  numbers with a mean of 1.0 and a standard deviation of 0.5

  Many libraries only generate uniformly distributed random numbers. If so, you may use 
  one of these algorithms 
  [https://en.wikipedia.org/wiki/Normal_distribution#Generating_values_from_normal_distribution]
  """

  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

import util.

main => go.


go =>
   _ = 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.
go => true.


% 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)).