Back to dada's perl lab
|
|
Statistical Moments |
Back to the Win32 Shootout Back to dada's perl lab |
| 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 $
// http://dada.perl.it/shootout/
//
// Transliterated from the Java implementation.
//
// contributed by Isaac Gouy
using System;
using System.Collections;
namespace LanguageShootout
{
class Moments
{
[STAThread]
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;
numbers.Add(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);
Readln(p^);
list.Add(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) ->
{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
Readln(input, num);
list.Add( Pointer(num) );
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
line max-line r@ read-line throw
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
\ others, addr1 u2 and addr3 u3 are the two partitions.
\ 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 }
addr u 1- 2/ floats + addr u 2/ floats +
addr addr u 1- floats + 2over quantile
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
-- read the file
main = do input <- getContents
putAns (lines input)
-- print out the answers
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
where conv x = fst (head (readFloat x))
-- 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 $
-- http://dada.perl.it/shootout/
-- Ada 95 code by C.C.
-- Annotated Ada Reference Manual ISO/IEC 8652:1995: http://www.ada-auth.org/
with System, Ada.Numerics.Generic_Elementary_Functions;
with Ada.Unchecked_Deallocation, Ada.Text_IO;
procedure Moments is
type Real is digits Positive'Max (15, System.Max_Digits);
package AF is new
Ada.Numerics.Generic_Elementary_Functions (Float_Type => Real);
package Io renames Ada.Text_IO;
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/
;;; from Brad Knotwell
(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)
(do ((line (read-line) (read-line)))
((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 $
// http://dada.perl.it/shootout
//
// 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 {
BufferedReader in = new BufferedReader(new InputStreamReader(System.in));
while ((line = in.readLine()) != null) {
num = Double.parseDouble(line);
sum += num;
nums.add(new Double(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
local line = read()
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 -->
io__read_file_as_string(_Res, Contents),
{ 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
import sumcol; // reuse bufferedStdInReader()
void main(String[] args){
double num, sum = 0;
ArrayList |
| 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
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
|
| 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
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
|
| 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()
{
array input = Stdio.stdin.read()/"\n";
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 $
# http://dada.perl.it/shootout/
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
line := STDIN readline
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)
(setf end (read-sequence buffer stream))
(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
nums = map(atof, sys.stdin.readlines())
sum = reduce(operator.add, nums)
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
pad = " "
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
pad = "0"
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
s = left(strip(s), width, pad)
else
s = right(strip(s), width, pad)
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
STDIN.readlines().each{|line|
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)
read_values
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]
-- Values read from stdin
read_values is
local
value: INTEGER
index: INTEGER
do
from
index := values.lower
std_input.read_line
until std_input.end_of_input loop
value := std_input.last_string.to_integer
values.force(value, index)
std_input.read_line
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 $
% http://dada.perl.it/shootout/
%
% 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
set nums [read stdin]
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
NumBlob = WScript.StdIn.ReadAll
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
// [dada] 2001-09-19 had to change power to pow for VC++
#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
Readln(input, num);
list.Add( Pointer(num) );
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.
|