stat: procedure options (main); /* 21 May 2014 */

stats: procedure (values, mean, standard_deviation);
   declare (values(*), mean, standard_deviation) float;
   declare n fixed binary (31) initial ( (hbound(values,1)) );

   mean = sum(values)/n;

   standard_deviation = sqrt( sum(values - mean)**2 / n);

end stats;

   declare values (*) float controlled;
   declare (mean, stddev) float;
   declare bin(0:9) fixed;
   declare (i, n) fixed binary (31);

   do n = 100, 1000, 10000, 100000;
      allocate values(n);
      values = random();
      call stats (values, mean, stddev);

      if n = 100 then
         do;
            bin = 0;
            do i = 1 to 100;
               bin(10*values(i)) += 1;
            end;
            put skip list ('Histogram for 100 values:');
            do i = 0 to 9;  /* display histogram */
               put skip list (repeat('.', bin(i)) );
            end;
         end;

      put skip list (n || ' values: mean=' || mean, 'stddev=' || stddev);
      free values;
   end;

end stat;