Statistical Moments Back to the Win32 Shootout Back to dada's perl lab

[The Original Shootout]   [NEWS]   [FAQ]   [Methodology]   [Platform Details]   [Acknowledgements]   [Scorecard]
All Source For Statistical Moments
moments.awka
```# \$Id: moments.gawk,v 1.1 2001/05/13 07:09:06 doug Exp \$
# http://www.bagley.org/~doug/shootout/

BEGIN {
delete ARGV;
sum = 0;
n = 0;
}

{
nums[n++] = \$1;
sum += \$1;
}

END {
mean = sum/n;
for (num in nums) {
dev = nums[num] - mean;
if (dev > 0) { avg_dev += dev; } else { avg_dev -= dev; }
vari += dev^2;
skew += dev^3;
kurt += dev^4;
}
avg_dev /= n;
vari /= (n - 1);
std_dev = sqrt(vari);

if (vari > 0) {
skew /= (n * vari * std_dev);
kurt = kurt/(n * vari * vari) - 3.0;
}

nums[n] = nums[0];
heapsort(n, nums);

mid = int(n/2)+1;
median = (n % 2) ? nums[mid] : (nums[mid] + nums[mid-1])/2;

printf("n:                  %d\n", n);
printf("median:             %f\n", median);
printf("mean:               %f\n", mean);
printf("average_deviation:  %f\n", avg_dev);
printf("standard_deviation: %f\n", std_dev);
printf("variance:           %f\n", vari);
printf("skew:               %f\n", skew);
printf("kurtosis:           %f\n", kurt);
}

function heapsort (n, ra) {
l = int(0.5+n/2) + 1
ir = n;
for (;;) {
if (l > 1) {
rra = ra[--l];
} else {
rra = ra[ir];
ra[ir] = ra[1];
if (--ir == 1) {
ra[1] = rra;
return;
}
}
i = l;
j = l * 2;
while (j <= ir) {
if (j < ir && ra[j] < ra[j+1]) { ++j; }
if (rra < ra[j]) {
ra[i] = ra[j];
j += (i = j);
} else {
j = ir + 1;
}
}
ra[i] = rra;
}
}
```
moments.bcc
```/* -*- mode: c -*-
* \$Id: moments.gcc,v 1.7 2001/01/05 02:07:51 doug Exp \$
* http://www.bagley.org/~doug/shootout/
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAXLINELEN 128

int
compare_doubles (const double *a, const double *b) {
return (int) (*a - *b);
}

int
main() {
char line[MAXLINELEN];
int i, n = 0, mid = 0;
double sum = 0.0;
double mean = 0.0;
double average_deviation = 0.0;
double standard_deviation = 0.0;
double variance = 0.0;
double skew = 0.0;
double kurtosis = 0.0;
double median = 0.0;
double deviation = 0.0;
int array_size = 4096;

double *nums = (double *)malloc(array_size * sizeof(double));

while (fgets(line, MAXLINELEN, stdin)) {
sum += (nums[n++] = atof(line));
if (n == array_size) {
array_size *= 2;
nums = (double *)realloc(nums, array_size * sizeof(double));
}
}
mean = sum/n;
for (i=0; i<n; i++) {
deviation = nums[i] - mean;
average_deviation += fabs(deviation);
variance += pow(deviation,2);
skew += pow(deviation,3);
kurtosis += pow(deviation,4);
}
average_deviation /= n;
variance /= (n - 1);
standard_deviation = sqrt(variance);
if (variance) {
skew /= (n * variance * standard_deviation);
kurtosis = (kurtosis/(n * variance * variance)) - 3.0;
}

qsort(nums, n, sizeof (double),
(int (*)(const void *, const void *))compare_doubles);
mid = (n/2);
median = n % 2 ? nums[mid] : (nums[mid] + nums[mid-1])/2;

free(nums);

printf("n:                  %d\n", n);
printf("median:             %f\n", median);
printf("mean:               %f\n", mean);
printf("average_deviation:  %f\n", average_deviation);
printf("standard_deviation: %f\n", standard_deviation);
printf("variance:           %f\n", variance);
printf("skew:               %f\n", skew);
printf("kurtosis:           %f\n", kurtosis);
return(0);
}
```
moments.csharp
```// \$Id: moments.csharp,v 1.0 2002/11/27 13:19:00 dada Exp \$
//
// Transliterated from the Java implementation.
//
// contributed by Isaac Gouy

using System;
using System.Collections;

namespace LanguageShootout
{
class Moments
{
static void Main(string[] args)
{
String line;
ArrayList numbers = new ArrayList();
double num, sum = 0.0;
double mean = 0.0;
double average_deviation = 0.0;
double standard_deviation = 0.0;
double variance = 0.0;
double skew = 0.0;
double kurtosis = 0.0;
double median = 0.0;
double deviation = 0.0;
int i, n, mid = 0;

while ( (line = Console.ReadLine()) != null ) {
num = Double.Parse(line);
sum += num;
}

n = numbers.Count;
mean = sum / n;
for (i=0; i<n; i++) {
deviation = (double)numbers[i] - mean;
average_deviation += Math.Abs(deviation);
variance += Math.Pow(deviation,2);
skew += Math.Pow(deviation,3);
kurtosis += Math.Pow(deviation,4);
}
average_deviation /= n;
variance /= (n - 1);
standard_deviation = Math.Sqrt(variance);
if (variance != 0.0) {
skew /= (n * variance * standard_deviation);
kurtosis = kurtosis/(n * variance * variance) - 3.0;
}

numbers.Sort();

mid = n / 2;
median = (n % 2 != 0) ?
(double)numbers[mid] :
((double)numbers[mid] + (double)numbers[mid-1]) / 2;

Console.WriteLine("n:                  {0:d}", n);
Console.WriteLine("median:             {0:f6}", median);
Console.WriteLine("mean:               {0:f6}", mean);
Console.WriteLine("average_deviation:  {0:f6}", average_deviation);
Console.WriteLine("standard_deviation: {0:f6}", standard_deviation);
Console.WriteLine("variance:           {0:f6}", variance);
Console.WriteLine("skew:               {0:f6}", skew);
Console.WriteLine("kurtosis:           {0:f6}", kurtosis);
}
}
}

```
moments.cygperl
```#!/usr/local/bin/perl
# \$Id: moments.perl,v 1.5 2001/01/05 22:36:44 doug Exp \$
# http://www.bagley.org/~doug/shootout/

use strict;

my @nums = <STDIN>;
my \$sum = 0;
foreach (@nums) { \$sum += \$_ }
my \$n = scalar(@nums);
my \$mean = \$sum/\$n;
my \$average_deviation = 0;
my \$standard_deviation = 0;
my \$variance = 0;
my \$skew = 0;
my \$kurtosis = 0;
foreach (@nums) {
my \$deviation = \$_ - \$mean;
\$average_deviation += abs(\$deviation);
\$variance += \$deviation**2;
\$skew += \$deviation**3;
\$kurtosis += \$deviation**4;
}
\$average_deviation /= \$n;
\$variance /= (\$n - 1);
\$standard_deviation = sqrt(\$variance);

if (\$variance) {
\$skew /= (\$n * \$variance * \$standard_deviation);
\$kurtosis = \$kurtosis/(\$n * \$variance * \$variance) - 3.0;
}

@nums = sort { \$a <=> \$b } @nums;
my \$mid = int(\$n/2);
my \$median = (\$n % 2) ? \$nums[\$mid] : (\$nums[\$mid] + \$nums[\$mid-1])/2;

printf("n:                  %d\n", \$n);
printf("median:             %f\n", \$median);
printf("mean:               %f\n", \$mean);
printf("average_deviation:  %f\n", \$average_deviation);
printf("standard_deviation: %f\n", \$standard_deviation);
printf("variance:           %f\n", \$variance);
printf("skew:               %f\n", \$skew);
printf("kurtosis:           %f\n", \$kurtosis);
```
moments.delphi
```program moments;

uses Classes, math;

function Compare(A, B : Pointer) : integer;
begin
if PDouble(A)^>PDouble(B)^ then
Result := 1
else if PDouble(A)^<PDouble(B)^ Then
Result:=-1
else
Result:=0;
end;

var
i, N, middle : integer;
list : TList;
median, mean, deviation,
average_deviation, standard_deviation,
variance, skew, kurtosis, sum: double;
p: PDouble;
begin
list := TList.Create;
While Not Eof(input) do begin
new(p);
end;
N := list.Count;
sum:=0.0;
for i:=0 to N-1 do
sum:=sum+PDouble(list.items[i])^;

mean := sum / N;
average_deviation:=0.0;
variance:=0.0;
skew:=0.0;
kurtosis:=0.0;

for i:=0 to N-1 do begin
deviation:=PDouble(list.items[i])^ - mean;
average_deviation:=average_deviation + Abs(deviation);
variance:=variance+power(deviation,2);
skew:=skew+power(deviation,3);
kurtosis:=kurtosis+power(deviation,4);
end;
average_deviation := average_deviation / N;
variance := variance / (N-1);
standard_deviation := Sqrt(variance);

if variance<>0 then begin
skew := skew / (N * variance * standard_deviation);
kurtosis := kurtosis / (N * variance * variance ) - 3.0;
end;

list.Sort(Compare);
middle := N Div 2;
if (N Mod 2) <> 0 then
median:=PDouble(list.items[middle])^
else
median:=(PDouble(list.items[middle])^+PDouble(list.items[middle-1])^)*0.5;

WriteLn('n:                  ', N);
WriteLn('median:             ', median:6:6);
WriteLn('mean:               ', mean:6:6);
WriteLn('average_deviation:  ', average_deviation:6:6);
WriteLn('standard_deviation: ', standard_deviation:6:6);
WriteLn('variance:           ', variance:6:6);
WriteLn('skew:               ', skew:6:6);
WriteLn('kurtosis:           ', kurtosis:6:6);
end.

```
moments.erlang
```%% \$Id: moments.erl,v 1.0 2002/09/24 12:09:00 dada Exp \$
%%
%% Code contributed by Isaac Gouy
%%
%% Usage: start from command line with
%%     erlc moments.erl
%%     erl -noinput -s moments main 200000

-module(moments).
-export([main/0]).
-import(lists, [sort/1, foldl/3, nth/2]).

accumulate([], Mean, Ad, Av, As, Ak) ->
accumulate([H|T], Mean, Ad, Av, As, Ak) ->
D = H - Mean,
D2 = D * D,
accumulate(T, Mean, Ad + abs(D), Av + D2, As + (D2 * D), Ak + (D2 * D2)).

median(L, N) -> medianS(sort(L), N).

medianS(L, N) ->
Mid = N div 2,
case N rem 2 of
0 -> (nth(Mid, L) + nth(Mid + 1, L)) / 2;
1 -> nth(Mid, L)
end.

skew(N, As, V, SD) when V > 0.0 -> As / (N * V * SD);
skew(N, As, V, SD) -> 0.

kurtosis(N, Ak, V) when V > 0.0 -> Ak / (N * V * V) - 3;
kurtosis(N, Ak, V) -> 0.

stats(L) ->
N = length(L),
io:fwrite("~s~w~n", ["n:                  ", N]),
io:fwrite("~s~f~n", ["median:             ", median(L, N)]),

Mean = foldl(fun(X, Sum) -> X + Sum end, 0.0, L) / N,
io:fwrite("~s~f~n", ["mean:               ", Mean]),

{Ad, Av, As, Ak} = accumulate(L, Mean, 0.0, 0.0, 0.0, 0.0),
io:fwrite("~s~f~n", ["average_deviation:  ", Ad / N]),

Variance = Av / (N - 1),
Standard_deviation = math:sqrt(Variance),
io:fwrite("~s~f~n", ["standard_deviation: ", Standard_deviation]),
io:fwrite("~s~f~n", ["variance:           ", Variance]),
io:fwrite("~s~f~n", ["skew:               ", skew(N, As, Variance, Standard_deviation)]),
io:fwrite("~s~f~n", ["kurtosis:           ", kurtosis(N, Ak, Variance)]),
io:fwrite("~n ", []).

main() ->
stats(lists:seq(1,500,1)),
halt(0).

```
moments.fpascal
```Program moments;
uses SysUtils, Classes;

function Power(Base : Real ; Exponent: Integer): Real;
var i : integer;
var pow : real;
begin
pow := Base;
For i:= 2 To Exponent do pow := pow * Base;
Power := pow;
end;

function Compare(A, B : Pointer) : longint;
begin
if A > B then
Compare := 1
else if A < B Then
Compare := -1
else
Compare := 0;
end;

var
i, N, sum, num, middle : longint;
list : TList;
median, mean, deviation,
average_deviation, standard_deviation,
variance, skew, kurtosis : real;
begin
list := TList.Create;
While Not Eof(input) do
begin
end;
N := list.Count;
For i := 0 To N-1 do Inc(sum, longint(list.Items[i]));
mean := sum / N;
average_deviation := 0;
standard_deviation := 0;
variance := 0;
skew := 0;
kurtosis := 0;

For i := 0 To N-1 do
begin
deviation := longint(list.Items[i]) - mean;
average_deviation := average_deviation + Abs(deviation);
variance := variance + Power(deviation, 2);
skew := skew + Power(deviation, 3);
kurtosis := kurtosis + Power(deviation, 4);

end;
average_deviation := average_deviation / N;
variance := variance / (N-1);
standard_deviation := Sqrt(variance);

If variance <> 0 Then
begin
skew := skew / (N * variance * standard_deviation);
kurtosis := kurtosis / (N * variance * variance ) - 3.0;
end;

list.Sort(@Compare);

middle := N Div 2;

If (N Mod 2) <> 0 Then
median := longint(list.Items[middle])
Else
median := (longint(list.Items[middle]) + longint(list.Items[middle-1])) / 2;

WriteLn('n:                  ', N);
WriteLn('median:             ', median:6:6);
WriteLn('mean:               ', mean:6:6);
WriteLn('average_deviation:  ', average_deviation:6:6);
WriteLn('standard_deviation: ', standard_deviation:6:6);
WriteLn('variance:           ', variance:6:6);
WriteLn('skew:               ', skew:6:6);
WriteLn('kurtosis:           ', kurtosis:6:6);
end.
```
moments.gawk
```# \$Id: moments.gawk,v 1.1 2001/05/13 07:09:06 doug Exp \$
# http://www.bagley.org/~doug/shootout/

BEGIN {
delete ARGV;
sum = 0;
n = 0;
}

{
nums[n++] = \$1;
sum += \$1;
}

END {
mean = sum/n;
for (num in nums) {
dev = nums[num] - mean;
if (dev > 0) { avg_dev += dev; } else { avg_dev -= dev; }
vari += dev^2;
skew += dev^3;
kurt += dev^4;
}
avg_dev /= n;
vari /= (n - 1);
std_dev = sqrt(vari);

if (vari > 0) {
skew /= (n * vari * std_dev);
kurt = kurt/(n * vari * vari) - 3.0;
}

nums[n] = nums[0];
heapsort(n, nums);

mid = int(n/2)+1;
median = (n % 2) ? nums[mid] : (nums[mid] + nums[mid-1])/2;

printf("n:                  %d\n", n);
printf("median:             %f\n", median);
printf("mean:               %f\n", mean);
printf("average_deviation:  %f\n", avg_dev);
printf("standard_deviation: %f\n", std_dev);
printf("variance:           %f\n", vari);
printf("skew:               %f\n", skew);
printf("kurtosis:           %f\n", kurt);
}

function heapsort (n, ra) {
l = int(0.5+n/2) + 1
ir = n;
for (;;) {
if (l > 1) {
rra = ra[--l];
} else {
rra = ra[ir];
ra[ir] = ra[1];
if (--ir == 1) {
ra[1] = rra;
return;
}
}
i = l;
j = l * 2;
while (j <= ir) {
if (j < ir && ra[j] < ra[j+1]) { ++j; }
if (rra < ra[j]) {
ra[i] = ra[j];
j += (i = j);
} else {
j = ir + 1;
}
}
ra[i] = rra;
}
}
```
moments.gcc
```/* -*- mode: c -*-
* \$Id: moments.gcc,v 1.7 2001/01/05 02:07:51 doug Exp \$
* http://www.bagley.org/~doug/shootout/
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAXLINELEN 128

int
compare_doubles (const double *a, const double *b) {
return (int) (*a - *b);
}

int
main() {
char line[MAXLINELEN];
int i, n = 0, mid = 0;
double sum = 0.0;
double mean = 0.0;
double average_deviation = 0.0;
double standard_deviation = 0.0;
double variance = 0.0;
double skew = 0.0;
double kurtosis = 0.0;
double median = 0.0;
double deviation = 0.0;
int array_size = 4096;

double *nums = (double *)malloc(array_size * sizeof(double));

while (fgets(line, MAXLINELEN, stdin)) {
sum += (nums[n++] = atof(line));
if (n == array_size) {
array_size *= 2;
nums = (double *)realloc(nums, array_size * sizeof(double));
}
}
mean = sum/n;
for (i=0; i<n; i++) {
deviation = nums[i] - mean;
average_deviation += fabs(deviation);
variance += pow(deviation,2);
skew += pow(deviation,3);
kurtosis += pow(deviation,4);
}
average_deviation /= n;
variance /= (n - 1);
standard_deviation = sqrt(variance);
if (variance) {
skew /= (n * variance * standard_deviation);
kurtosis = (kurtosis/(n * variance * variance)) - 3.0;
}

qsort(nums, n, sizeof (double),
(int (*)(const void *, const void *))compare_doubles);
mid = (n/2);
median = n % 2 ? nums[mid] : (nums[mid] + nums[mid-1])/2;

free(nums);

printf("n:                  %d\n", n);
printf("median:             %f\n", median);
printf("mean:               %f\n", mean);
printf("average_deviation:  %f\n", average_deviation);
printf("standard_deviation: %f\n", standard_deviation);
printf("variance:           %f\n", variance);
printf("skew:               %f\n", skew);
printf("kurtosis:           %f\n", kurtosis);
return(0);
}
```
moments.gforth
```\ -*- mode: forth -*-
\ \$Id: moments.gforth,v 1.1 2001/06/03 12:10:25 doug Exp \$
\ http://www.bagley.org/~doug/shootout/
\ from Anton Ertl

1024 constant max-line
create line max-line 2 + allot

: input-floats
>r 0e begin
while
line swap >float 0= abort" float expected"
fdup f, f+
repeat
rdrop drop ;

: compute-loop
dup 0 d>f fdup { f: n } f/ { f: mean }
0e fdup fdup fdup
floats bounds do {  f: avg-deviation f: variance f: skew f: kurtosis }
i f@ mean f- { f: deviation }
deviation fabs avg-deviation f+
deviation fdup f* fdup variance f+
fswap deviation f* fdup skew f+
fswap deviation f* kurtosis f+
float +loop
frot n 1e f- f/ to variance
frot to avg-deviation
variance fsqrt { f: standard-deviation }
variance f0<> if
n variance fdup f* f* f/ 3e f-
fswap n variance f*  standard-deviation f* f/ fswap
endif
fswap variance standard-deviation avg-deviation n f/ mean ;

float- -1 floats ,

: partition { first last -- last-smaller first-larger }
\ partition array addr1 u1 into all elements less than pivot and all
\ lessthan-xt ( elemptr1 elemptr2 -- f ) compares the two elements
first last + 1 rshift faligned f@ { f: pivot }
first last begin
begin
pivot dup f@ f< over first u> and
while
float-
repeat
swap begin
dup last u< over f@ pivot f< and
while
float+
repeat
2dup u>=
while
dup f@ over f@ dup f! over f!
float+ swap float-
repeat ;

: quantile  recursive
\ sorts the array [first,last] such that the contained part of
\ [quant-low,quant-high] is the same as in the fully sorted array.
{ quant-low quant-high }
begin { first last }
first quant-high u< quant-low last u< and
while
first last partition
last quant-low quant-high quantile
first swap
repeat ;

: median { addr u -- rmedian }
f@ f@ f+ f2/ ;

: ff.
f\$ dup >r 0<
IF '0 emit ELSE scratch r@ min type r@ precision - zeros ENDIF
'. emit
r@ negate zeros
scratch r> 0 max /string 0 max type ;

create nums \ s" moments.input" r/o open-file throw input-floats
stdin input-floats
nums here over - float /
." n:                  " dup 0 .r cr
compute-loop
nums here over - float / median  9 set-precision
." median:             " ff. cr
." mean:               " ff. cr
." average_deviation:  " ff. cr
." standard_deviation: " ff. cr 11 set-precision
." variance:           " ff. cr  7 set-precision
." skew:               " ff. cr
." kurtosis:           " ff. cr
bye
```
moments.ghc
```-- \$Id: moments.ghc,v 1.1 2001/02/14 22:12:21 doug Exp \$
-- http://www.bagley.org/~doug/shootout/
-- from Brian Gregor

module Main where

import IO
import System
import Numeric

main = do input <- getContents
putAns (lines input)

putAns :: [String] -> IO ()
putAns st_nums = do
putStrLn ("n:                  " ++ (showInt (truncate n) ""))
putStrLn ("median:             " ++ (showFFloat (Just 6) (median nums n) ""))
putStrLn ("mean:               " ++ (showFFloat (Just 6) mean ""))
putStrLn ("average_deviation:  " ++ (showFFloat (Just 6) avg_dev ""))
putStrLn ("standard_deviation: " ++ (showFFloat (Just 6) std_dev ""))
putStrLn ("variance:           " ++ (showFFloat (Just 6) var ""))
putStrLn ("skew:               " ++ (showFFloat (Just 6) skew ""))
putStrLn ("kurtosis:           " ++ (showFFloat (Just 6) kurt ""))
where
n = fromIntegral (length nums)
nums = strToDoub st_nums
mean = (sum nums) / n
deviation = [x-mean | x <- nums]
avg_dev = (sum [abs x | x <- deviation])/ n
var = (sum [x**2 | x <- deviation]) / (n-1)
std_dev = sqrt var
skew = (sum [x**3 | x <- deviation]) / (n*var*std_dev)
kurt = (sum [x**4 | x <- deviation]) / (n*var*var)-3.0

-- convert the strings to doubles
strToDoub :: [String] -> [Double]
strToDoub nums = map conv nums

-- calculate the median
median :: [Double] -> Double -> Double
median nums n = mid (mSort nums)
where
mid x
| odd (length x) = x!! midpt
| otherwise       = ((x!!(midpt-1)) + (x!!midpt)) / 2.0
midpt :: Int
midpt = floor (n/2)

-- Sorting: the various languages use various algorithms
-- here's  an optimized mergesort from
-- "Algorithms - a Functional Approach" by
-- Fethi Rabhe & Guy Lapalme
split :: (Ord a) => [a] -> [[a]]
split [] = []
split (x:xs) = [x]:split xs

merge :: (Ord a) => [a] -> [a] -> [a]
merge [] b  = b
merge a [] = a
merge a@(x:xs) b@(y:ys)
| (x<=y) = x : (merge xs b)
| otherwise = y : (merge a ys)

mergepairs :: (Ord a) => [[a]] -> [[a]]
mergepairs [] = []
mergepairs x@[l] = x
mergepairs (l1:l2:rest) = (merge l1 l2) : (mergepairs \$! rest)

-- The actual sort
mSort :: (Ord a) => [a] -> [a]
mSort l = ms (split l)
where  ms [r] = r
ms l   = ms (mergepairs l)
```
moments.gnat
```-- \$Id: moments.gnat,v 1.0 2003/06/11 12:08:00 dada Exp \$
-- Ada 95 code by C.C.

procedure Moments is
type Real is digits Positive'Max (15, System.Max_Digits);
package AF is new
package Rio is new Ada.Text_IO.Float_IO (Num => Real);
procedure Put (Item : Real; Fore : Io.Field := 0; Aft : Io.Field := 6;
Exp : Io.Field := 0) renames Rio.Put;

generic
type Item_Type is private;      --  Component type of array to be sorted
with function "<=" (X, Y : Item_Type) return Boolean;
type Sequence is array (Integer range <>) of Item_Type;
package Sort_Pkg is
procedure Quick_Sort (S : in out Sequence);
end Sort_Pkg;  --  Copied from Southampton Ada MPICH supercomputer bindings

package body Sort_Pkg is
procedure Quick_Sort (S : in out Sequence) is
procedure Sort_Recursive (Lwb, Upb : Integer) is
Pivot    : Item_Type := S (Upb);
Front    : Integer := Lwb;
Back     : Integer := Upb;
Temp     : Item_Type;
begin
if Lwb < Upb then
while  (Front <= Back)  loop
while not  (Pivot <= S (Front)) loop
Front := Front + 1;
end loop;
while not  (S (Back) <= Pivot) loop
Back := Back - 1;
end loop;
if Front <= Back then
Temp := S (Front);
S (Front) := S (Back);
S (Back) := Temp;
Front := Front + 1;
Back  := Back - 1;
end if;
end loop;
Sort_Recursive (Lwb, Back);
Sort_Recursive (Front, Upb);
end if;
end Sort_Recursive;
begin
Sort_Recursive (S'First, S'Last);
end Quick_Sort;
end Sort_Pkg;

type Real_Array is array (Integer range <>) of Real;
type Real_Array_Ptr is access Real_Array;
procedure Deallocate_Real_Array is new Ada.Unchecked_Deallocation (
Object => Real_Array, Name => Real_Array_Ptr);
package Sort is new Sort_Pkg (Real, "<=" => "<=", Sequence => Real_Array);

Data           : Real_Array_Ptr := new Real_Array (1 .. 0);
Temp_Array     : Real_Array_Ptr;
Dev, D_2, Mean, Median     : Real;
Standard_Deviation         : Real;
Sum, Avg_Abs_Deviation     : Real := 0.0;
Variance, Skew, Kurtosis   : Real := 0.0;
M              : Natural := 0;
begin
while not Io.End_Of_File loop
M := M + 1;
if M > Data'Last then                  --  Lengthen array of strings
Temp_Array := new Real_Array (
1 .. 4 + Positive (Real'Ceiling (1.70 * Real (M))));
Temp_Array (Data'Range) := Data.all;
Deallocate_Real_Array (Data);
Data := Temp_Array;     --  Replace the old array with the new one
end if;
Rio.Get (Item => Data (M));      -- Read  numbers between whitespace
Sum := Sum + Data (M);
end loop;
Mean := Sum / Real (M);
for K in 1 .. M loop
Dev := Data (K) - Mean;
Avg_Abs_Deviation := Avg_Abs_Deviation + abs Dev;
D_2 := Dev * Dev;
Variance := Variance + D_2;
Skew := Skew + (D_2 * Dev);
Kurtosis := Kurtosis + (D_2 * D_2);
end loop;
Avg_Abs_Deviation := Avg_Abs_Deviation / Real (M);
Variance := Variance / Real (M - 1);
Standard_Deviation := AF.Sqrt (Variance);
if Variance < 10.0 * Real'Model_Epsilon then
Io.Put_Line ("> Reduced accuracy results: 0 = ((Variance/10 + 1) - 1)");
else
Skew := Skew / (Real (M) * Variance * Standard_Deviation);
Kurtosis := -3.0 + Kurtosis / (Real (M) * Variance * Variance);
end if;
Sort.Quick_Sort (S => Data (1 .. M));
if 1 = (M mod 2) then
Median := Data ((M + 1) / 2);
else
Median := (Data (M / 2) + Data (1 + M / 2)) / 2.0;
end if;
Io.Put_Line ("n:                 " & Integer'Image (M));
Io.Put ("median:             "); Put (Median);              Io.New_Line;
Io.Put ("mean:               "); Put (Mean);                Io.New_Line;
Io.Put ("average_deviation:  "); Put (Avg_Abs_Deviation);   Io.New_Line;
Io.Put ("standard_deviation: "); Put (Standard_Deviation);  Io.New_Line;
Io.Put ("variance:           "); Put (Variance);            Io.New_Line;
Io.Put ("skew:               "); Put (Skew);                Io.New_Line;
Io.Put ("kurtosis:           "); Put (Kurtosis);            Io.New_Line;
end Moments;

```
moments.guile
```#!/usr/local/bin/guile \
-e main -s
!#

;;; \$Id: moments.guile,v 1.2 2001/06/29 23:12:37 doug Exp \$
;;; http://www.bagley.org/~doug/shootout/

(use-modules (ice-9 format))

(define sum 0)
(define nums '())
(define (compute-all mean n continuation)
(let ((average-deviation 0) (standard-deviation 0) (variance 0) (skew 0)
(kurtosis 0) (mid 0) (median 0) (deviation 0) (tmp-lst nums))
(do ((num (car tmp-lst) (if (eq? tmp-lst '()) '() (car tmp-lst))))
((eq? num '()) (begin (set! variance (/ variance (1- n)))
(set! standard-deviation (sqrt variance))
(if (>; variance 0.0)
(begin (set! skew (/ skew (* n variance standard-deviation)))
(set! kurtosis (- (/ kurtosis (* n variance variance)) 3))))
(set! nums (sort-list nums >;))
(set! mid (/ n 2))
(set! median (if (= (remainder n 2) 0)
(/ (+ (list-ref nums mid)
(list-ref nums (1- mid)))
2)
(list-ref nums mid)))
(continuation n median mean
(/ average-deviation n)
standard-deviation variance
skew kurtosis)))
(let ((deviation (- num mean)))
(begin (set! tmp-lst (cdr tmp-lst))
(set! average-deviation (+ average-deviation (abs deviation)))
(set! variance (+ variance (expt deviation 2)))
(set! skew (+ skew (expt deviation 3)))
(set! kurtosis (+ kurtosis (expt deviation 4))))))))

(define output-format
"n:                  ~D
median:             ~,6F
mean:               ~,6F
average_deviation:  ~,6F
standard_deviation: ~,6F
variance:           ~,6F
skew:               ~,6F
kurtosis:           ~,6F
")

(define (main args)
((eof-object? line)
(compute-all (/ sum (length nums)) (length nums)
(lambda (x . y) (display (apply format (cons output-format (cons x y)))))))
(let ((num (string->;number line)))
(begin (set! nums (cons num nums)) (set! sum (+ sum num))))))
```
moments.ici
```// \$Id: moments.ici,v 1.0 2003/01/03 11:55:00 dada Exp \$
//
// contributed by Tim Long

sum := 0.0;
nums := array();
forall (f in gettokens(stdin, "\n", ""))
{
push(nums, f = float(f));
sum += f;
}

n := nels(nums);
mean := sum / n;

deviation := 0.0;
average_deviation := 0.0;
standard_deviation := 0.0;
variance := 0.0;
skew := 0.0;
kurtosis := 0.0;

forall (num in nums)
{
deviation = num - mean;
average_deviation += abs(deviation);
variance += (t := deviation * deviation);
skew += (t *= deviation);
kurtosis += (t *= deviation);
}
average_deviation /= n;
variance /= (n - 1);
standard_deviation = sqrt(variance);

if (variance > 0.0)
{
skew /= n * variance * standard_deviation;
kurtosis = kurtosis / (n * variance * variance) - 3.0;
}

sort(nums);
mid := n / 2;
if (n % 2 == 0)
median = (nums[mid] + nums[mid - 1])/2;
else
median = nums[mid];

printf("n:                  %d\n", n);
printf("median:             %f\n", median);
printf("mean:               %f\n", mean);
printf("average_deviation:  %f\n", average_deviation);
printf("standard_deviation: %f\n", standard_deviation);
printf("variance:           %f\n", variance);
printf("skew:               %f\n", skew);
printf("kurtosis:           %f\n", kurtosis);
```
moments.java
```// \$Id: moments.java,v 1.2 2001/01/05 01:35:56 doug Exp \$
// http://www.bagley.org/~doug/shootout/

import java.io.*;
import java.util.*;
import java.text.*;
import java.lang.Math;

public class moments {
public static void main(String[] args) {
String line;
Vector nums = new Vector();
double num, sum = 0.0;
double mean = 0.0;
double average_deviation = 0.0;
double standard_deviation = 0.0;
double variance = 0.0;
double skew = 0.0;
double kurtosis = 0.0;
double median = 0.0;
double deviation = 0.0;
int i, n, mid = 0;

try {
while ((line = in.readLine()) != null) {
num = Double.parseDouble(line);
sum += num;
}
} catch (IOException e) {
System.err.println(e);
return;
}

n = nums.size();
mean = sum/n;
for (i=0; i<n; i++) {
deviation = ((Double)nums.get(i)).doubleValue() - mean;
average_deviation += Math.abs(deviation);
variance += Math.pow(deviation,2);
skew += Math.pow(deviation,3);
kurtosis += Math.pow(deviation,4);
}
average_deviation /= n;
variance /= (n - 1);
standard_deviation = Math.sqrt(variance);
if (variance != 0.0) {
skew /= (n * variance * standard_deviation);
kurtosis = kurtosis/(n * variance * variance) - 3.0;
}

Collections.sort(nums);

mid = (n/2);
median = (n % 2 != 0) ?
((Double)nums.get(mid)).doubleValue() :
(((Double)nums.get(mid)).doubleValue() +
((Double)nums.get(mid-1)).doubleValue())/2;

NumberFormat nf = NumberFormat.getInstance();
nf.setMaximumFractionDigits(13);
nf.setGroupingUsed(false);
nf.setMaximumFractionDigits(6);
nf.setMinimumFractionDigits(6);

System.out.println("n:                  " + n);
System.out.println("median:             " + nf.format(median));
System.out.println("mean:               " + nf.format(mean));
System.out.println("average_deviation:  " + nf.format(average_deviation));
System.out.println("standard_deviation: " + nf.format(standard_deviation));
System.out.println("variance:           " + nf.format(variance));
System.out.println("skew:               " + nf.format(skew));
System.out.println("kurtosis:           " + nf.format(kurtosis));
}
}
```
moments.lcc
```/* -*- mode: c -*-
* \$Id: moments.gcc,v 1.7 2001/01/05 02:07:51 doug Exp \$
* http://www.bagley.org/~doug/shootout/
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAXLINELEN 128

int
compare_doubles (const double *a, const double *b) {
return (int) (*a - *b);
}

int
main() {
char line[MAXLINELEN];
int i, n = 0, mid = 0;
double sum = 0.0;
double mean = 0.0;
double average_deviation = 0.0;
double standard_deviation = 0.0;
double variance = 0.0;
double skew = 0.0;
double kurtosis = 0.0;
double median = 0.0;
double deviation = 0.0;
int array_size = 4096;

double *nums = (double *)malloc(array_size * sizeof(double));

while (fgets(line, MAXLINELEN, stdin)) {
sum += (nums[n++] = atof(line));
if (n == array_size) {
array_size *= 2;
nums = (double *)realloc(nums, array_size * sizeof(double));
}
}
mean = sum/n;
for (i=0; i<n; i++) {
deviation = nums[i] - mean;
average_deviation += fabs(deviation);
variance += pow(deviation,2);
skew += pow(deviation,3);
kurtosis += pow(deviation,4);
}
average_deviation /= n;
variance /= (n - 1);
standard_deviation = sqrt(variance);
if (variance) {
skew /= (n * variance * standard_deviation);
kurtosis = (kurtosis/(n * variance * variance)) - 3.0;
}

qsort(nums, n, sizeof (double),
(int (*)(const void *, const void *))compare_doubles);
mid = (n/2);
median = n % 2 ? nums[mid] : (nums[mid] + nums[mid-1])/2;

free(nums);

printf("n:                  %d\n", n);
printf("median:             %f\n", median);
printf("mean:               %f\n", mean);
printf("average_deviation:  %f\n", average_deviation);
printf("standard_deviation: %f\n", standard_deviation);
printf("variance:           %f\n", variance);
printf("skew:               %f\n", skew);
printf("kurtosis:           %f\n", kurtosis);
return(0);
}
```
moments.lua
```-- \$Id: moments.lua,v 1.2 2001/01/05 01:35:56 doug Exp \$
-- http://www.bagley.org/~doug/shootout/
-- implemented by: Roberto Ierusalimschy

local nums = {}
local n = 0
local sum = 0
while 1 do
if line == nil then break end
line = line+0        -- convert line to number
sum = sum + line
n = n + 1
nums[n] = line
end

local mean = sum/n

local average_deviation, variance, skew, kurtosis = 0, 0, 0, 0

for i = 1, n do
local deviation = nums[i] - mean
average_deviation = average_deviation + abs(deviation)
variance = variance + deviation^2
skew = skew + deviation^3
kurtosis = kurtosis + deviation^4
end

average_deviation = average_deviation/n
variance = variance/(n-1)
local standard_deviation = sqrt(variance)
if variance ~= 0 then
skew = skew / (n * variance * standard_deviation)
kurtosis = kurtosis/(n * variance * variance) - 3.0
end

sort(nums)
local mid = floor(n/2)
local median
if mod(n,2) == 1 then
median = nums[mid+1]
else
median = (nums[mid] + nums[mid+1])/2
end

write(format("n:                  %d\n", n))
write(format("median:             %f\n", median))
write(format("mean:               %f\n", mean))
write(format("average_deviation:  %f\n", average_deviation))
write(format("standard_deviation: %f\n", standard_deviation))
write(format("variance:           %f\n", variance))
write(format("skew:               %f\n", skew))
write(format("kurtosis:           %f\n", kurtosis))
```
moments.lua5
```-- \$Id: moments.lua,v 1.2 2001/01/05 01:35:56 doug Exp \$
-- http://www.bagley.org/~doug/shootout/
-- contributed by Roberto Ierusalimschy

local nums = {}
local n = 0
local sum = 0
for line in io.lines() do
line = line+0        -- convert line to number
sum = sum + line
n = n + 1
nums[n] = line
end

local mean = sum/n

local average_deviation, variance, skew, kurtosis = 0, 0, 0, 0

for i = 1, n do
local deviation = nums[i] - mean
average_deviation = average_deviation + math.abs(deviation)
variance = variance + deviation^2
skew = skew + deviation^3
kurtosis = kurtosis + deviation^4
end

average_deviation = average_deviation/n
variance = variance/(n-1)
local standard_deviation = math.sqrt(variance)
if variance ~= 0 then
skew = skew / (n * variance * standard_deviation)
kurtosis = kurtosis/(n * variance * variance) - 3.0
end

table.sort(nums)
local mid = math.floor(n/2)
local median
if math.mod(n,2) == 1 then
median = nums[mid+1]
else
median = (nums[mid] + nums[mid+1])/2
end

io.write(string.format("n:                  %d\n", n))
io.write(string.format("median:             %f\n", median))
io.write(string.format("mean:               %f\n", mean))
io.write(string.format("average_deviation:  %f\n", average_deviation))
io.write(string.format("standard_deviation: %f\n", standard_deviation))
io.write(string.format("variance:           %f\n", variance))
io.write(string.format("skew:               %f\n", skew))
io.write(string.format("kurtosis:           %f\n", kurtosis))

```
moments.mawk
```# \$Id: moments.mawk,v 1.1 2001/05/13 07:09:06 doug Exp \$
# http://www.bagley.org/~doug/shootout/

BEGIN {
delete ARGV;
sum = 0;
n = 0;
}

{
nums[n++] = \$1;
sum += \$1;
}

END {
mean = sum/n;
for (num in nums) {
dev = nums[num] - mean;
if (dev > 0) { avg_dev += dev; } else { avg_dev -= dev; }
vari += dev^2;
skew += dev^3;
kurt += dev^4;
}
avg_dev /= n;
vari /= (n - 1);
std_dev = sqrt(vari);

if (vari > 0) {
skew /= (n * vari * std_dev);
kurt = kurt/(n * vari * vari) - 3.0;
}

nums[n] = nums[0];
heapsort(n, nums);

mid = int(n/2)+1;
median = (n % 2) ? nums[mid] : (nums[mid] + nums[mid-1])/2;

printf("n:                  %d\n", n);
printf("median:             %f\n", median);
printf("mean:               %f\n", mean);
printf("average_deviation:  %f\n", avg_dev);
printf("standard_deviation: %f\n", std_dev);
printf("variance:           %f\n", vari);
printf("skew:               %f\n", skew);
printf("kurtosis:           %f\n", kurt);
}

function heapsort (n, ra) {
l = int(0.5+n/2) + 1
ir = n;
for (;;) {
if (l > 1) {
rra = ra[--l];
} else {
rra = ra[ir];
ra[ir] = ra[1];
if (--ir == 1) {
ra[1] = rra;
return;
}
}
i = l;
j = l * 2;
while (j <= ir) {
if (j < ir && ra[j] < ra[j+1]) { ++j; }
if (rra < ra[j]) {
ra[i] = ra[j];
j += (i = j);
} else {
j = ir + 1;
}
}
ra[i] = rra;
}
}
```
moments.mercury
```%% \$Id: moments.mercury,v 1.1 2001/07/29 15:25:39 doug Exp \$
%% http://www.bagley.org/~doug/shootout/
%% from Fergus Henderson

:- module mytest.
:- interface.
:- import_module io.

:- pred main(io__state, io__state).
:- mode main(di, uo) is det.

:- implementation.
:- import_module array, string, float, math, int, list, require.

main -->
{ Lines = string__words((pred('\n'::in) is semidet), Contents) },
{ Count = length(Lines) },
{ array__init(Count, 0.0, Array0) },
{ count_and_sum(Lines, 0, 0.0, Array0, _Count2, Sum, Array) },
{ Mean = Sum / float(Count) },
process(0, Count, Mean, 0.0, 0.0, 0.0, 0.0, Array).

:- pred count_and_sum(list(string), int, float, array(float),
int, float, array(float)).
:- mode count_and_sum(in, in, in, array_di, out, out, array_uo) is det.
count_and_sum([], Count, Sum, Array, Count, Sum, Array).
count_and_sum([L|Ls], Count0, Sum0, Array0, Count, Sum, Array) :-
(if string__to_float(L, V) then Val = V else error("float conversion")),
count_and_sum(Ls, Count0 + 1, Sum0 + Val, Array0^elem(Count0) := Val,
Count, Sum, Array).

:- pred process(int, int, float, float, float, float, float, array(float),
io__state, io__state).
:- mode process(in, in, in, in, in, in, in, array_di, di, uo) is det.
process(I, Count, Mean,
SumAbsDeviations, SumVariance, SumSkew, SumKurtosis, Array0) -->
(if { I < Count } then
{ Val = Array0 ^ elem(I) },
{ Dev = Val - Mean },
{ Dev2 = Dev * Dev },
{ Dev3 = Dev2 * Dev },
{ Dev4 = Dev2 * Dev2 },
process(I + 1, Count, Mean, SumAbsDeviations + abs(Dev),
SumVariance + Dev2, SumSkew + Dev3,
SumKurtosis + Dev4, Array0)
else
{
AverageDeviation = SumAbsDeviations / float(Count),
Variance = SumVariance / float(Count - 1),
StandardDeviation = sqrt(Variance),
(if Variance \= 0.0 then
Skew = SumSkew / (float(Count) * Variance *
StandardDeviation),
Kurtosis = (SumKurtosis / (float(Count) *
Variance * Variance)) - 3.0
else
Skew = 0.0,
Kurtosis = 0.0
),
Array = sort(Array0),
Mid = (Count//2),
Median = (if Count rem 2 = 1 then Array^elem(Mid)
else (Array^elem(Mid) + Array^elem(Mid - 1)) / 2.0)
},
format("n:                  %d\n", [i(Count)]),
format("median:             %f\n", [f(Median)]),
format("mean:               %f\n", [f(Mean)]),
format("average_deviation:  %f\n", [f(AverageDeviation)]),
format("standard_deviation: %f\n", [f(StandardDeviation)]),
format("variance:           %f\n", [f(Variance)]),
format("skew:               %f\n", [f(Skew)]),
format("kurtosis:           %f\n", [f(Kurtosis)])
).
```
moments.mingw32
```/* -*- mode: c -*-
* \$Id: moments.gcc,v 1.7 2001/01/05 02:07:51 doug Exp \$
* http://www.bagley.org/~doug/shootout/
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAXLINELEN 128

int
compare_doubles (const double *a, const double *b) {
return (int) (*a - *b);
}

int
main() {
char line[MAXLINELEN];
int i, n = 0, mid = 0;
double sum = 0.0;
double mean = 0.0;
double average_deviation = 0.0;
double standard_deviation = 0.0;
double variance = 0.0;
double skew = 0.0;
double kurtosis = 0.0;
double median = 0.0;
double deviation = 0.0;
int array_size = 4096;

double *nums = (double *)malloc(array_size * sizeof(double));

while (fgets(line, MAXLINELEN, stdin)) {
sum += (nums[n++] = atof(line));
if (n == array_size) {
array_size *= 2;
nums = (double *)realloc(nums, array_size * sizeof(double));
}
}
mean = sum/n;
for (i=0; i<n; i++) {
deviation = nums[i] - mean;
average_deviation += fabs(deviation);
variance += pow(deviation,2);
skew += pow(deviation,3);
kurtosis += pow(deviation,4);
}
average_deviation /= n;
variance /= (n - 1);
standard_deviation = sqrt(variance);
if (variance) {
skew /= (n * variance * standard_deviation);
kurtosis = (kurtosis/(n * variance * variance)) - 3.0;
}

qsort(nums, n, sizeof (double),
(int (*)(const void *, const void *))compare_doubles);
mid = (n/2);
median = n % 2 ? nums[mid] : (nums[mid] + nums[mid-1])/2;

free(nums);

printf("n:                  %d\n", n);
printf("median:             %f\n", median);
printf("mean:               %f\n", mean);
printf("average_deviation:  %f\n", average_deviation);
printf("standard_deviation: %f\n", standard_deviation);
printf("variance:           %f\n", variance);
printf("skew:               %f\n", skew);
printf("kurtosis:           %f\n", kurtosis);
return(0);
}
```
moments.nice
```/* The Great Win32 Language Shootout http://dada.perl.it/shootout/
contributed by Isaac Gouy (Nice novice)

To compile:
nicec --sourcepath=.. -d=. -a moments.jar moments

To run:
java -jar moments.jar < input.txt > out.txt
*/

import java.io.*;
import java.text.*;

import random; // reuse floatFormat

void main(String[] args){
double num, sum = 0;
ArrayList numbers = new ArrayList();
try {
?String line;
num = Double.parseDouble(line);
sum += num;
};
}
catch (IOException e) {
System.err.println(e);
}

double mean = 0.0;
double average_deviation = 0.0;
double standard_deviation = 0.0;
double variance = 0.0;
double skew = 0.0;
double kurtosis = 0.0;
double median = 0.0;
double deviation = 0.0;
int i, n, mid = 0;

n = numbers.size();
mean = sum/n;
for (i=0; i
```
moments.ocaml
```(*
* \$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
!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
```
moments.ocamlb
```(*
* \$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
!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
```
moments.perl
```#!/usr/local/bin/perl
# \$Id: moments.perl,v 1.5 2001/01/05 22:36:44 doug Exp \$
# http://www.bagley.org/~doug/shootout/

use strict;

my @nums = <STDIN>;
my \$sum = 0;
foreach (@nums) { \$sum += \$_ }
my \$n = scalar(@nums);
my \$mean = \$sum/\$n;
my \$average_deviation = 0;
my \$standard_deviation = 0;
my \$variance = 0;
my \$skew = 0;
my \$kurtosis = 0;
foreach (@nums) {
my \$deviation = \$_ - \$mean;
\$average_deviation += abs(\$deviation);
\$variance += \$deviation**2;
\$skew += \$deviation**3;
\$kurtosis += \$deviation**4;
}
\$average_deviation /= \$n;
\$variance /= (\$n - 1);
\$standard_deviation = sqrt(\$variance);

if (\$variance) {
\$skew /= (\$n * \$variance * \$standard_deviation);
\$kurtosis = \$kurtosis/(\$n * \$variance * \$variance) - 3.0;
}

@nums = sort { \$a <=> \$b } @nums;
my \$mid = int(\$n/2);
my \$median = (\$n % 2) ? \$nums[\$mid] : (\$nums[\$mid] + \$nums[\$mid-1])/2;

printf("n:                  %d\n", \$n);
printf("median:             %f\n", \$median);
printf("mean:               %f\n", \$mean);
printf("average_deviation:  %f\n", \$average_deviation);
printf("standard_deviation: %f\n", \$standard_deviation);
printf("variance:           %f\n", \$variance);
printf("skew:               %f\n", \$skew);
printf("kurtosis:           %f\n", \$kurtosis);
```
moments.pike
```#!/usr/local/bin/pike// -*- mode: pike -*-
// \$Id: moments.pike,v 1.6 2001/01/05 01:35:56 doug Exp \$
// http://www.bagley.org/~doug/shootout/
// from: Fredrik Noring

class Moments
{
int N;
float median;
float mean;
float average_deviation;
float standard_deviation;
float variance;
float skew;
float kurtosis;

void create(array(float) v)
{
float sum = `+(@v);
N = sizeof(v);
mean = sum / N;

foreach(v, float i)
{
float deviation = i - mean;
average_deviation += abs(deviation);
variance += pow(deviation, 2);
skew += pow(deviation, 3);
kurtosis += pow(deviation, 4);
}

average_deviation /= N;
variance /= (N - 1);
standard_deviation = sqrt(variance);

if (variance)
{
skew /= (N * variance * standard_deviation);
kurtosis = kurtosis/(N * variance * variance) - 3.0;
}

sort(v);
int mid = N/2;
median = N % 2 ? v[mid] : (v[mid] + v[mid-1])/2;
}
};

int main()
{
Moments m=Moments( (array(float)) input[..sizeof(input)-2] );

write("n:                  %d\n", m->N);
write("median:             %.6f\n", m->median);
write("mean:               %.6f\n", m->mean);
write("average_deviation:  %.6f\n", m->average_deviation);
write("standard_deviation: %.6f\n", m->standard_deviation);
write("variance:           %.6f\n", m->variance);
write("skew:               %.6f\n", m->skew);
write("kurtosis:           %.6f\n", m->kurtosis);
}
```
moments.pliant
```# \$Id: moments.pliant,v 1.0 2002/02/08 12:43:00 dada Exp \$

module "/pliant/language/context.pli"
module "/pliant/language/stream.pli"
module "/pliant/language/stream/pipe.pli"

function heapsort n ra
arg Int n ; arg_rw Array:Int ra
var Int rra
var Int i
var Int j
var Int l := (n\2) + 1
var Int ir := n

part heapsort_loop
if l>1
l := l - 1
rra := ra:l
else
rra := ra:ir
ra:ir := ra:1
ir := ir - 1
if ir=1
ra:1 := rra
leave heapsort_loop
i := l
j := l*2

while j<=ir
if j<ir and ra:j < ra:(j+1)
j := j + 1
if rra < ra:j
ra:i := ra:j
i := j
j := j + i
else
j := ir + 1
ra:i := rra
restart heapsort_loop

gvar Str line := ""
(gvar Stream STDIN) open "handle:0" in
gvar Array:Int nums
gvar uInt sum := 0
gvar Int i
gvar Int n
gvar Int mid
gvar Float mean
gvar Float deviation
gvar Float average_deviation := 0
gvar Float standard_deviation := 0
gvar Float variance := 0
gvar Float devpow3
gvar Float skew := 0
gvar Float kurtosis := 0
gvar Float median

while (STDIN atend) = false
i := 0
if(line parse i any)
nums += i

for i 0 nums:size-1
sum += nums:i

n := nums:size
mean := sum/n

for i 0 nums:size-1
deviation := nums:i - mean
average_deviation := average_deviation + (abs deviation)
variance := variance + ((abs deviation) ^ 2.0)
devpow3 := (abs deviation) ^ 3.0
if deviation < 0
devpow3 := -devpow3
skew := skew + devpow3
kurtosis := kurtosis + ((abs deviation) ^ 4.0)

average_deviation := average_deviation / n
variance := variance / (n-1)
standard_deviation := variance^0.5

if variance <> 0
skew := skew / (n * variance * standard_deviation)
kurtosis := kurtosis / (n * variance * variance) - 3.0

heapsort n nums
mid := n\2
median := shunt (n%2>0) nums:mid (nums:mid+nums:mid-1)/2

console "n:                  " n eol
console "median:             " (string median "fixed 6") eol
console "mean:               " (string mean "fixed 6") eol
console "average_deviation:  " (string average_deviation "fixed 6") eol
console "standard_deviation: " (string standard_deviation "fixed 6") eol
console "variance:           " (string variance "fixed 6") eol
console "skew:               " (string skew "fixed 6") eol
console "kurtosis:           " (string kurtosis "fixed 6") eol
```
moments.poplisp
```;;; -*- mode: lisp -*-
;;; \$Id: moments.poplisp,v 1.0 2002/05/08 17:37:00 dada Exp \$

(declaim (optimize (speed 3) (debug 0) (safety 0) (space 0) (compilation-speed 0)))

(defun quicksort (vec lo hi) ;; modified from from Roger Corman's posting in cll
(declare (fixnum hi lo) (type (simple-array double-float) vec))
(if (>; hi lo)
(let* ((mid (round (+ lo hi) 2))
(i lo)
(j (+ hi 1))
(p (aref vec mid)))
(declare (fixnum i j) (double-float p))
(rotatef (aref vec mid) (aref vec lo)) ;; swap mid element to first
(loop
(loop do (incf i)
until (or (>; i hi) (> p (aref vec i))))
(loop do (decf j)
until (or (<;= j lo) (> (aref vec j) p)))
(if (<; j i) (return))
(rotatef (aref vec i)(aref vec j)))
(rotatef (aref vec lo) (aref vec j)) ;;put partition element in place
(quicksort vec lo (- j 1))  (quicksort vec i hi))) vec)

(defun do-moments (data n mean)
(declare (fixnum n) (double-float mean) (type (simple-array double-float) data))
(let ((average_deviation 0.0d0)
(standard_deviation 0.0d0)
(variance 0.0d0)
(skew 0.0d0)
(kurtosis 0.0d0)
(median 0.0d0))
(declare (double-float mean average_deviation standard_deviation
variance skew kurtosis median))
(declare (inline quicksort))
(loop for i fixnum from 0 below n do
(let* ((deviation (- (the double-float (aref data i)) mean))
(deviation2 (* deviation deviation))
(deviation3 (* deviation deviation2))
(deviation4 (* deviation deviation3)))
(incf average_deviation (abs deviation))
(incf variance deviation2)
(incf skew deviation3)
(incf kurtosis deviation4)))

(setq average_deviation (/ average_deviation n))
(setq variance (/ variance (1- n)))
(setq standard_deviation (sqrt variance))

(cond ((>; variance 0.0d0)
(setq skew (/ skew (* n variance standard_deviation)))
(setq kurtosis (- (/ kurtosis (* (coerce n 'double-float)
variance variance))
3.0d0))))
(setf data (quicksort data 0 (1- n)))
(let ((mid (/ n 2)))
(fixnum mid)
(if (zerop (mod n 2))
(setq median (/ (+ (the double-float (aref data mid))
(the double-float (aref data (1- mid))))
2.0d0))
(setq median (aref data mid))))
(format t "n:                  ~A~%" n)
(format t "median:             ~,6K~%" median)
(format t "mean:               ~,6K~%" mean)
(format t "average_deviation:  ~,6K~%" average_deviation)
(format t "standard_deviation: ~,6K~%" standard_deviation)
(format t "variance:           ~,6K~%" variance)
(format t "skew:               ~,6K~%" skew)
(format t "kurtosis:           ~,6K~%" kurtosis)))

(let ((buffer (make-string 4096))
(start 0)
(end 0)
(result 0.0d0)
(char #\X)
(stream *standard-input*)
(eof-p nil))
(declare (fixnum start end) (double-float result))
(labels ((get-char ()
(when (= start end)
(setf start 0)
(when (zerop end)
(setf eof-p t)
(setf char #\Z) ;any non-digit will do
(return-from get-char nil)))
(setf char (schar buffer start))
(incf start))
(get-dfloat ();; parse double float hack someone should rewrite this
(let ((minusp nil)
(expminusp nil)
(before-dp 0)
(after-dp 0)
(dec-digits 0)
(exponent 0))
(declare (fixnum before-dp after-dp exponent dec-digits)
(inline digit-char-p char=))
(loop while (and
(not
(or (and (char= #\- char)
(setq minusp t))
(digit-char-p char 10)))
(get-char)))
(loop
do (let ((weight (digit-char-p char 10)))
(declare (type (or null fixnum) weight))
(if weight
(setq before-dp (+ weight (the fixnum (* before-dp 10))))
(return)))
until (not (get-char)))
(if minusp (setf before-dp (- before-dp)))
(when (char= #\. char)
(loop while (get-char)
do (let ((weight (digit-char-p char 10)))
(declare (type (or null (signed-byte 32)) weight))
(if weight
(setq after-dp (+ weight (the fixnum (* after-dp 10)))
dec-digits (the fixnum (1+ dec-digits)))
(return)))))
(when (or (char= #\e char) (char= #\E char))
(get-char)
(when (char= #\- char)
(setq expminusp t)
(get-char))
(loop
do (let ((weight (digit-char-p char 10)))
(declare (type (or null fixnum) weight))
(if weight
(setq exponent (+ weight (the fixnum (* exponent 10))))
(return)))
until (not (get-char)))
(if expminusp (setf exponent (- exponent))))
(setq result
(float (*
(+ (float before-dp 1.0d0)
(if (zerop after-dp) 0.0d0
(* (float after-dp 1.0d0)
(if (zerop dec-digits) 1.0d0
(expt 10.0d0 (float (- dec-digits) 1.0d0))))))
(if (zerop exponent) 1.0d0
(expt 10.0d0 (float exponent 1.0d0)))) 1.0d0)))))

(let ((sum 0.0d0)
nums )
(declare (double-float sum) (inline vector-push-extend))
(let* ((array-size 10000)
(numbuffer (make-array array-size :element-type 'double-float))
(buflist (list numbuffer)) ;; Doug's idea put these together later
(fill-pointer 0))
(loop
(get-dfloat)
(if (not eof-p)
(progn
(incf sum result)
(setf (aref numbuffer fill-pointer) result)
(incf fill-pointer)
(when (= fill-pointer array-size)
(push
(setf numbuffer (make-array array-size :element-type 'double-float))
buflist)
(setf fill-pointer 0)))
(return)))
(let* ((num-arrays (length buflist))
(num-elem (+ (* (1- num-arrays) array-size) fill-pointer)))
(setf nums (make-array  num-elem :element-type 'double-float))
(locally (declare (type (simple-array double-float) nums))
(loop for i fixnum from 0 to (1- num-arrays) do
(setf (subseq nums (* i array-size))
(the (simple-array double-float)
(elt buflist (- (1- num-arrays) i))))) ;;buflist is rev'd
(do-moments nums num-elem (/ sum num-elem))))))))
```
moments.python
```#!/usr/local/bin/python
# \$Id: moments.python,v 1.4 2001/05/09 00:58:40 doug Exp \$
# http://www.bagley.org/~doug/shootout/

import sys, string, math, operator

def main():
sum = 0
nums = []
atof = string.atof

n = len(nums)
mean = sum/n
average_deviation = 0
standard_deviation = 0
variance = 0
skew = 0
kurtosis = 0

for num in nums:
deviation = num - mean
average_deviation += abs(deviation)
variance += deviation**2;
skew += deviation**3;
kurtosis += deviation**4

average_deviation /= n
variance /= (n - 1)
standard_deviation = math.sqrt(variance)

if variance > 0.0:
skew /= (n * variance * standard_deviation)
kurtosis = kurtosis/(n * variance * variance) - 3.0

nums.sort()
mid = n / 2

if (n % 2) == 0:
median = (nums[mid] + nums[mid-1])/2
else:
median = nums[mid]

print "n:                  %d" % n
print "median:             %f" % median
print "mean:               %f" % mean
print "average_deviation:  %f" % average_deviation
print "standard_deviation: %f" % standard_deviation
print "variance:           %f" % variance
print "skew:               %f" % skew
print "kurtosis:           %f" % kurtosis

main()

```
moments.rexx
```NUMERIC DIGITS 10

n = 1
nums. = 0
sum = 0
DO WHILE LINES() <> 0
PARSE LINEIN L
IF L <> "" THEN DO
INTERPRET "sum = sum + " L
INTERPRET "nums." || n || "=" || L
n = n + 1
END
END

n = n - 1
nums.0 = n
mean = sum/n
average_deviation = 0
standard_deviation = 0
variance = 0
skew = 0
kurtosis = 0
DO I = 1 TO nums.0
deviation = nums.I - mean
average_deviation = average_deviation + ABS(deviation)
variance = variance + deviation ** 2
skew = skew + deviation ** 3
kurtosis = kurtosis + deviation ** 4
END
average_deviation = average_deviation / n
variance = variance / (n-1)
standard_deviation = SQRT(variance)

IF variance <> 0 THEN DO
skew = skew / (n * variance * standard_deviation)
kurtosis = kurtosis / (n * variance * variance) - 3
END

call qqsort 1, n

mid = n%2+1
if n//2 <> 0 then do
median = nums.mid
end
else do
pmid = mid - 1
median = (nums.mid + nums.pmid) / 2
end

SAY sprintf("n:                  %d", n)
SAY sprintf("median:             %-.6f", median)
SAY sprintf("mean:               %-.6f", mean)
SAY sprintf("average_deviation:  %-.6f", average_deviation)
SAY sprintf("standard_deviation: %-.6f", standard_deviation)
SAY sprintf("variance:           %-.6f", variance)
SAY sprintf("skew:               %-.6f", skew)
SAY sprintf("kurtosis:           %-.6f", kurtosis)

EXIT

SQRT: PROCEDURE
ARG n

ans = n / 2
prevans = n
do until prevans = ans
prevans = ans
ans = ( prevans + ( n / prevans ) ) / 2
end
return ans

SQRTOLD:
PROCEDURE
PARSE ARG n
parse value format(n,,,,0) with mant "E" exp
if exp = "" then exp = 0
if exp//2 < > 0 then do
mant = mant * 10
exp = exp -1
end
root = 0
do 10
do digit = 9 by -1 to 0,
while,
(root + digit) ** 2 > mant
end
root = root + digit
if root**2 = mant then leave
root = root * 10
mant = mant * 100
exp = exp -2
end
return root * 10**(exp/2)

sprintf: procedure
argno = 1
string = ""
start = 1
len = length(arg(1))

do until(p >= len)
s = ""
argno = argno + 1
p = pos("%", arg(1), start)
if p = 0 then
do
p = len + 1
end
if substr(arg(1), p, 1) == "%" then
do
s = "%"
end
string = string || substr(arg(1), start, p - start)
start = p + 1
p = verify(arg(1), "%cdfsx", "M", start)
if p = 0 then
leave
spec = substr(arg(1), start, p - start + 1)
start = p + 1
r = right(spec, 1)
spec = delstr(spec, length(spec), 1)
if left(spec,1) == "-" then
do
left = 1
spec = substr(spec, 2)
end
else
do
left = 0
spec = substr(spec, 1)
end
if spec \== "" then
parse var spec width "." prec
else
do
width = 0
prec = 0
end
if \datatype(width, "W") then
width = 0
if \datatype(prec, "W") then
prec = 0

select

when r == "s" then
do
if width = 0 then
width = length(arg(argno))
if prec \= 0 then
s = left(arg(argno), prec)
else
s = arg(argno)
end

when r == "d" then
do
if width = 0 then
width = length(arg(argno))
s = format(arg(argno), length(arg(argno)), 0)
end

when r == "f" then
do
if arg(argno) > -1 & arg(argno) < 1 then
parse value arg(argno) with int "." frac
if width = 0 & prec = 0 then
do
d = 1
if arg(argno) < 0 then d = 2
width = digits() + d
prec = digits() - (length(int)) + d - 1
end
if width = 0 then
width = len - prec
s = format(arg(argno), width, prec, 0)
end

when r == "x" then
do
if width = 0 then
width = length(arg(argno))
s = d2x(arg(argno))
if prec \= 0 then
s = left(s, prec)
end

when r == "%" then
do
argno = argno - 1
end

otherwise
nop

end

if r \== "%" then
do
if left then
else
end
string = string || s
end
return string

Fast Quick sort

/*
This is a fast quick sort routine.

Author: Ruediger Wilke
*/

qqsort: procedure expose nums.

arg lf, re

if re -lf < 9 then
do lf = lf to re -1

m = lf

do j = lf +1 to re
if nums.j < nums.m then
m = j
end

t = nums.m; nums.m = nums.lf; nums.lf = t

end
else
do
i = lf
j = re
k = (lf + re)%2
t = nums.k

do until i > j

do while nums.i < t
i = i + 1
end

do while nums.j > t
j = j - 1
end

if i <= j then
do
xchg = nums.i
nums.i = nums.j
nums.j = xchg
i = i + 1
j = j - 1
end

end

call qqsort lf, j
call qqsort i, re
end

return

```
moments.ruby
```#!/usr/local/bin/ruby
# -*- mode: ruby -*-
# \$Id: moments.ruby,v 1.5 2001/01/05 01:35:56 doug Exp \$
# http://www.bagley.org/~doug/shootout/

# throw away unused parameter sent by benchmark framework
ARGV.shift()

sum = 0.0
nums = []
num = nil
deviation = nil

num = Float(line)
nums << num
sum += num
}
n = nums.length()
mean = sum/n;
average_deviation = 0
standard_deviation = 0
variance = 0
skew = 0
kurtosis = 0

nums.each{|num|
deviation = num - mean
average_deviation += deviation.abs()
variance += deviation**2;
skew += deviation**3;
kurtosis += deviation**4
}
average_deviation /= n
variance /= (n - 1)
standard_deviation = Math.sqrt(variance)

if (variance > 0.0)
skew /= (n * variance * standard_deviation)
kurtosis = kurtosis/(n * variance * variance) - 3.0
end

nums.sort()
mid = n / 2

if (n % 2) == 0
median = (nums[mid] + nums[mid-1])/2
else
median = nums[mid]
end

printf("n:                  %d\n", n)
printf("median:             %f\n", median)
printf("mean:               %f\n", mean)
printf("average_deviation:  %f\n", average_deviation)
printf("standard_deviation: %f\n", standard_deviation)
printf("variance:           %f\n", variance)
printf("skew:               %f\n", skew)
printf("kurtosis:           %f\n", kurtosis)
```
moments.se
```-- -*- mode: eiffel -*-
-- \$Id: moments.se,v 1.2 2001/05/23 18:29:33 doug Exp \$
-- http://www.bagley.org/~doug/shootout/
-- from Steve Thompson

-- <LOC-OFF>
indexing
description: "This class performs the statistical moments test"
author : Steve Thompson
email  : "Steve_Thompson@prodigy.net"
date   : February 18, 2001
compile: "compile -clean -boost -no_split -O3 main.e -o main"
run    : "main < Input"
-- <LOC-ON>

class MOMENTS

creation make

feature -- Creation

make is
local
index: INTEGER
number: INTEGER
mid: INTEGER
deviation: DOUBLE
sorter: COLLECTION_SORTER[INTEGER]
do
!!values.make(0,499)
mean := sum / values.count

from index := values.lower until index > values.upper loop
number := values @ index
deviation := number - mean
average_deviation := average_deviation + deviation.abs
variance := variance + deviation ^ 2
skew := skew + deviation ^ 3
kurtosis := kurtosis + deviation ^ 4
index := index + 1
end

average_deviation := average_deviation / values.count
variance := variance / (values.count - 1)
standard_deviation := variance.sqrt --math.sqrt(variance)
if variance > 0.0 then
skew := skew / (values.count * variance * standard_deviation)
kurtosis := kurtosis / (values.count * variance * variance) - 3.0
end
sorter.sort(values)
mid := values.count // 2
if (values.count \\ 2) = 0 then
median := (values.item(mid) + values.item(mid - 1)) / 2
else
median := values @ mid
end

print ("n:                  " + values.count.to_string + "%N")
print ("median:             " + median.to_string + "%N")
print ("mean:               " + mean.to_string + "%N")
print ("average_deviation:  " + average_deviation.to_string + "%N")
print ("standard_deviation: " + standard_deviation.to_string + "%N")
print ("variance:           " + variance.to_string + "%N")
print ("skew:               " + skew.to_string + "%N")
print ("kurtosis:           " + kurtosis.to_string + "%N")
end -- make

feature -- Queries

sum: INTEGER

mean,
median,
average_deviation,
standard_deviation,
variance,
skew,
kurtosis: DOUBLE

values: ARRAY[INTEGER]

local
value: INTEGER
index: INTEGER
do

from
index := values.lower
until std_input.end_of_input loop
value := std_input.last_string.to_integer
values.force(value, index)
sum := sum + value
index := index + 1
end

end -- get_stdin

end
```
moments.slang
```% \$Id: moments.slang,v 1.0 2003/01/03 13:38:00 dada Exp \$
%
% contributed by John E. Davis

define main ()

{
variable nums = array_map (Double_Type, &atof, fgetslines (stdin));
variable n = length (nums);

variable sum = 0;
foreach (nums) sum += ();
variable mean = sum/n;

variable average_deviation = 0;
variable standard_deviation = 0;
variable variance = 0;
variable skew = 0;
variable kurtosis = 0;

% No slang programmer would code explicit loops.
foreach (nums)
{
variable num = ();
variable deviation = num - mean;
average_deviation += abs(deviation);
variance += deviation^2;
skew += deviation^3;
kurtosis += deviation^4;
}
average_deviation /= n;
variance /= (n - 1);
standard_deviation = sqrt(variance);

if (variance > 0.0)
{
skew /= (n * variance * standard_deviation);
kurtosis = kurtosis/(n * variance * variance) - 3.0;
}

nums = nums[array_sort(nums)];
variable mid = n/2;
variable median;

if (n mod 2)
median = nums[mid];
else
median = (nums[mid] + nums[mid-1])/2;

vmessage ("n:                  %d", n);
vmessage ("median:             %f", median);
vmessage ("mean:               %f", mean);
vmessage ("average_deviation:  %f", average_deviation);
vmessage ("standard_deviation: %f", standard_deviation);
vmessage ("variance:           %f", variance);
vmessage ("skew:               %f", skew);
vmessage ("kurtosis:           %f", kurtosis);
}
main ();
```
moments.smlnj
```(* -*- mode: sml -*-
* \$Id: moments.smlnj,v 1.3 2001/07/10 12:55:23 doug Exp \$
* http://www.bagley.org/~doug/shootout/
* from Stephen Weeks
* with help from Daniel Wang:
*   Modified to be more functional and use SML/NJ library sorting function
*)

structure Test : sig
val main : (string * string list) -> OS.Process.status
end = struct
val ins = TextIO.stdIn
fun loop (nums,sum) =
case TextIO.inputLine ins of
"" => (nums,sum)
| l => (case Real.fromString l of
NONE => raise Fail "invalid input"
| SOME num => loop(num::nums,sum+num))

fun printl [] = print "\n" | printl(h::t) = ( print h ; printl t )

fun r2s (x: real): string =
if Real.== (x, 0.0) then "0.000000"
else String.translate
(fn #"~" => "-" | c => str c)
(Real.fmt (StringCvt.FIX (SOME 6)) x)

fun main(name, args) =  let

val (nums,sum) = loop ([],0.0)
val nums = Array.fromList nums
val n = Array.length nums
val n_float = real n
val mean = sum / n_float

fun moments (x,{average_deviation,variance,skew,kurtosis}) = let
val deviation = x - mean
val average_deviation =
average_deviation + abs(deviation)
val dev2 = deviation * deviation
val variance = variance + dev2
val dev3 = dev2 * deviation
val skew = skew + dev3
val dev4 = dev3 * deviation
val kurtosis = kurtosis + dev4
in {average_deviation=average_deviation,
variance=variance,
skew=skew,
kurtosis=kurtosis}
end
val init = {average_deviation=0.0,
variance=0.0,
skew=0.0,
kurtosis=0.0}

val {average_deviation,variance,skew,kurtosis} =
Array.foldl moments init nums

val average_deviation = average_deviation / n_float
val variance = variance /  real (n - 1);
val standard_deviation = Real.Math.sqrt (variance)
val {skew,kurtosis} =
if variance > 0.0
then {skew=skew / n_float / variance / standard_deviation,
kurtosis=kurtosis / n_float / variance / variance - 3.0}
else {skew=skew,kurtosis=kurtosis}

val _ = ArrayQSort.sort Real.compare nums
val mid = Int.quot (n, 2)
val median =
if Int.rem (n, 2) = 1
then Array.sub (nums, mid)
else (Array.sub (nums, mid) +
Array.sub (nums, mid - 1)) / 2.0
in
printl ["n:                  ", Int.toString n, "\n",
"median:             ", r2s median, "\n",
"mean:               ", r2s mean, "\n",
"average_deviation:  ", r2s average_deviation, "\n",
"standard_deviation: ", r2s standard_deviation, "\n",
"variance:           ", r2s variance, "\n",
"skew:               ", r2s skew, "\n",
"kurtosis:           ", r2s kurtosis];
OS.Process.success
end
end

val _ = SMLofNJ.exportFn("moments", Test.main);
```
moments.tcl
```#!/usr/local/bin/tclsh
# \$Id: moments.tcl,v 1.3 2001/01/05 22:11:27 doug Exp \$
# http://www.bagley.org/~doug/shootout/

proc main {} {
set sum 0.0
foreach num \$nums {
set sum [expr {\$sum + \$num}]
}
set n [llength \$nums]
set mean [expr {\$sum / \$n}]
set average_deviation 0.0
set standard_deviation 0.0
set variance 0.0
set skew 0.0
set kurtosis 0.0

foreach num \$nums {
set deviation [expr {\$num - \$mean}]
set average_deviation [expr {\$average_deviation + abs(\$deviation)}]
set variance [expr {\$variance + pow(\$deviation, 2)}]
set skew [expr {\$skew + pow(\$skew, 3)}]
set kurtosis [expr {\$kurtosis + pow(\$deviation, 4)}]
}

set average_deviation [expr {\$average_deviation / \$n}]
set variance [expr {\$variance / (\$n - 1)}]
set standard_deviation [expr {sqrt(\$variance)}]

if {\$variance} {
set skew [expr {\$skew / (\$n * \$variance * \$standard_deviation)}]
set kurtosis [expr {\$kurtosis / (\$n * \$variance * \$variance) - 3.0}]
}

set nums [lsort -integer \$nums]
set mid [expr {int(\$n / 2)}]
if [expr {\$n % 2}] {
set median [lindex \$nums \$mid]
} else {
set a [lindex \$nums \$mid]
set b [lindex \$nums [expr {\$mid - 1}]]
set median [expr {(\$a + \$b) / 2.0}]
}

puts [format "n:                  %d" \$n]
puts [format "median:             %f" \$median]
puts [format "mean:               %f" \$mean]
puts [format "average_deviation:  %f" \$average_deviation]
puts [format "standard_deviation: %f" \$standard_deviation]
puts [format "variance:           %f" \$variance]
puts [format "skew:               %f" \$skew]
puts [format "kurtosis:           %f" \$kurtosis]
}

main
```
moments.vbscript
```<job>

<script language=JScript runat=server>

function SortNumeric(a, b) {
return ((+a > +b) ? 1 : ((+a < +b) ? -1 : 0));
}

function SortVBArray(arrVBArray) {
return arrVBArray.toArray().sort(SortNumeric).join('x');
}
</script>

<script language=VBScript>

Function SortArray(arrInput)
SortArray = Split(SortVBArray(arrInput), "x")
End Function

Sum = 0
N = 0
Num = Split(NumBlob, vbCrLf)
N = UBound(Num)

For A = 0 To N - 1
Sum = Sum + Num(A)
Next

mean = sum / N
average_deviation = 0
standard_deviation = 0
variance = 0
skew = 0
kurtosis = 0
For A = 0 To N - 1
deviation = Num(A) - mean
average_deviation = average_deviation + Abs(deviation)
variance = variance + deviation^2
skew = skew + deviation^3
kurtosis = kurtosis + deviation^4
Next
average_deviation = average_deviation / N
variance = variance / (N-1)
standard_deviation = Sqr(variance)
If variance Then
skew = skew / (N * variance * standard_deviation)
kurtosis = kurtosis / (n * variance * variance ) - 3.0
End If

SortNum = SortArray(Num)

middle = N/2 + 1

If (N Mod 2) Then
median = CInt(SortNum(middle))
Else
median = (CInt(SortNum(middle)) + CInt(SortNum(middle-1))) / 2
End If

WScript.Echo "n:                  " & N
WScript.Echo "median:             " & FormatNumber(median, 6, -1, 0, 0)
WScript.Echo "mean:               " & FormatNumber(mean, 6, -1, 0, 0)
WScript.Echo "average_deviation:  " & FormatNumber(average_deviation, 6, -1, 0, 0)
WScript.Echo "standard_deviation: " & FormatNumber(standard_deviation, 6, -1, 0, 0)
WScript.Echo "variance:           " & FormatNumber(variance, 6, -1, 0, 0)
WScript.Echo "skew:               " & FormatNumber(skew, 6, -1, 0, 0)
WScript.Echo "kurtosis:           " & FormatNumber(kurtosis, 6, -1, 0, 0)
</script>
</job>
```
moments.vc
```/* -*- mode: c -*-
* \$Id: moments.gcc,v 1.7 2001/01/05 02:07:51 doug Exp \$
* http://www.bagley.org/~doug/shootout/
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAXLINELEN 128

int
compare_doubles (const double *a, const double *b) {
return (int) (*a - *b);
}

int
main() {
char line[MAXLINELEN];
int i, n = 0, mid = 0;
double sum = 0.0;
double mean = 0.0;
double average_deviation = 0.0;
double standard_deviation = 0.0;
double variance = 0.0;
double skew = 0.0;
double kurtosis = 0.0;
double median = 0.0;
double deviation = 0.0;
int array_size = 4096;

double *nums = (double *)malloc(array_size * sizeof(double));

while (fgets(line, MAXLINELEN, stdin)) {
sum += (nums[n++] = atof(line));
if (n == array_size) {
array_size *= 2;
nums = (double *)realloc(nums, array_size * sizeof(double));
}
}
mean = sum/n;
for (i=0; i<n; i++) {
deviation = nums[i] - mean;
average_deviation += fabs(deviation);
variance += pow(deviation,2);
skew += pow(deviation,3);
kurtosis += pow(deviation,4);
}
average_deviation /= n;
variance /= (n - 1);
standard_deviation = sqrt(variance);
if (variance) {
skew /= (n * variance * standard_deviation);
kurtosis = (kurtosis/(n * variance * variance)) - 3.0;
}

qsort(nums, n, sizeof (double),
(int (*)(const void *, const void *))compare_doubles);
mid = (n/2);
median = n % 2 ? nums[mid] : (nums[mid] + nums[mid-1])/2;

free(nums);

printf("n:                  %d\n", n);
printf("median:             %f\n", median);
printf("mean:               %f\n", mean);
printf("average_deviation:  %f\n", average_deviation);
printf("standard_deviation: %f\n", standard_deviation);
printf("variance:           %f\n", variance);
printf("skew:               %f\n", skew);
printf("kurtosis:           %f\n", kurtosis);
return(0);
}
```
moments.vc++
```// -*- mode: c++ -*-
// \$Id
// http://www.bagley.org/~doug/shootout/
// Calculate statistical moments of a region, from Bill Lear

#include <iostream>
#include <vector>
#include <numeric>
#include <iterator>
#include <algorithm>
#include <cstdio>
#include <math.h>

using namespace std;

template <class T>
struct moments {
public:
template <class InputIterator>
moments(InputIterator begin, InputIterator end)
: median(0.0), mean(0.0), average_deviation(0.0),
standard_deviation(0.0), variance(0.0),
skew(0.0), kurtosis(0.0)
{
T sum = accumulate(begin, end, 0.0);
size_t N = end - begin;
mean = sum / N;

for (InputIterator i = begin; i != end; ++i) {
T deviation = *i - mean;
average_deviation += fabs(deviation);
variance += pow(deviation, 2);
skew += pow(deviation, 3);
kurtosis += pow(deviation, 4);
}

average_deviation /= N;
variance /= (N - 1);
standard_deviation = sqrt(variance);

if (variance) {
skew /= (N * variance * standard_deviation);
kurtosis = kurtosis/(N * variance * variance) - 3.0;
}

sort(begin, end);
size_t mid = N/2;
median = N % 2 ? *(begin+mid) : (*(begin+mid) + *(begin+mid-1))/2;
}

T median;
T mean;
T average_deviation;
T standard_deviation;
T variance;
T skew;
T kurtosis;
};

int main() {
vector<double> v;
char line[8];

ios_base::sync_with_stdio(false);
cin.tie(0);

while (cin.getline(line, 8)) {
v.push_back(atof(line));
}

moments<double> m(v.begin(), v.end());

int n = v.end() - v.begin();

printf("n:                  %d\n", n);
printf("median:             %f\n", m.median);
printf("mean:               %f\n", m.mean);
printf("average_deviation:  %f\n", m.average_deviation);
printf("standard_deviation: %f\n", m.standard_deviation);
printf("variance:           %f\n", m.variance);
printf("skew:               %f\n", m.skew);
printf("kurtosis:           %f\n", m.kurtosis);
}
```
moments.vpascal
```Program moments;

uses SysUtils, Classes;

function Power(Base : Real ; Exponent: Integer): Real;
var i : integer;
var pow : real;
begin
pow := Base;
For i:= 2 To Exponent do pow := pow * Base;
Power := pow;
end;

function Compare(A, B : Pointer) : integer;
begin
if Integer(A) > Integer(B) then
Result := 1
else if Integer(A) < Integer(B) Then
Result := -1
else
Result := 0;
end;

var
i, N, sum, num, middle : integer;
list : TList;
median, mean, deviation,
average_deviation, standard_deviation,
variance, skew, kurtosis : real;

begin
list := TList.Create;
While Not Eof(input) do
begin
end;
N := list.Count;
For i := 0 To N-1 do Inc(sum, Integer(list.Items[i]));
mean := sum / N;
average_deviation := 0;
standard_deviation := 0;
variance := 0;
skew := 0;
kurtosis := 0;

For i := 0 To N-1 do
begin
deviation := Integer(list.Items[i]) - mean;
average_deviation := average_deviation + Abs(deviation);
variance := variance + Power(deviation, 2);
skew := skew + Power(deviation, 3);
kurtosis := kurtosis + Power(deviation, 4);

end;
average_deviation := average_deviation / N;
variance := variance / (N-1);
standard_deviation := Sqrt(variance);

If variance <> 0 Then
begin
skew := skew / (N * variance * standard_deviation);
kurtosis := kurtosis / (N * variance * variance ) - 3.0;
end;

list.Sort(Compare);
middle := N Div 2;

If (N Mod 2) <> 0 Then
median := Integer(list.Items[middle])
Else
median := (Integer(list.Items[middle]) + Integer(list.Items[middle-1])) / 2;

WriteLn('n:                  ', N);
WriteLn('median:             ', median:6:6);
WriteLn('mean:               ', mean:6:6);
WriteLn('average_deviation:  ', average_deviation:6:6);
WriteLn('standard_deviation: ', standard_deviation:6:6);
WriteLn('variance:           ', variance:6:6);
WriteLn('skew:               ', skew:6:6);
WriteLn('kurtosis:           ', kurtosis:6:6);
end.
```