(*
 * $Id: moments.ocaml,v 1.9 2001/05/20 16:43:13 doug Exp $
 * http://www.bagley.org/~doug/shootout/
 * with help from Markus Mottl
 *)

let _ =
  let n = ref 0
  and num = ref 0.0
  and sum = ref 0.0
  and mean = ref 0.0
  and average_deviation = ref 0.0
  and standard_deviation = ref 0.0
  and variance = ref 0.0
  and skew = ref 0.0
  and kurtosis = ref 0.0
  and deviation = ref 0.0
  and size = ref 4096 in
  let nums_in = ref (Array.create !size 0.0) in

  try
    while true do
      num := read_float ();
      !nums_in.(!n) <- !num;
      sum := !sum +. !num;
      incr n;
      if !n = !size then begin
    nums_in := Array.append !nums_in (Array.create !size 0.0);
    size := !size * 2
      end
    done
  with End_of_file -> ();

  let nums = Array.create !n 0.0 in
  Array.blit !nums_in 0 nums 0 !n;

  let n_float = float_of_int !n in
  mean := !sum /. n_float;

  for i = 0 to !n - 1 do
    deviation := nums.(i) -. !mean;
    average_deviation := !average_deviation +. abs_float !deviation;
    let dev2 = !deviation *. !deviation in
    variance := !variance +. dev2;
    let dev3 = dev2 *. !deviation in
    skew := !skew +. dev3;
    let dev4 = dev3 *. !deviation in
    kurtosis := !kurtosis +. dev4;
  done;

  average_deviation := !average_deviation /. n_float;
  variance := !variance /. float_of_int (!n - 1);
  standard_deviation := sqrt !variance;

  if !variance > 0.0 then begin
    skew := !skew /. n_float /. !variance /. !standard_deviation;
    kurtosis := !kurtosis /. n_float /. !variance /. !variance -. 3.0;
  end;

  Array.stable_sort compare nums;

  let mid = !n lsr 1 in

  let median =
    if !n mod 2 = 1 then nums.(mid)
    else (nums.(mid) +. nums.(mid - 1)) /. 2.0 in

  Printf.printf "n:                  %d\n" !n;
  Printf.printf "median:             %f\n" median;
  Printf.printf "mean:               %f\n" !mean;
  Printf.printf "average_deviation:  %f\n" !average_deviation;
  Printf.printf "standard_deviation: %f\n" !standard_deviation;
  Printf.printf "variance:           %f\n" !variance;
  Printf.printf "skew:               %f\n" !skew;
  Printf.printf "kurtosis:           %f\n" !kurtosis