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