1%%%  This code was developped by IDEALX (http://IDEALX.org/) and
2%%%  contributors (their names can be found in the CONTRIBUTORS file).
3%%%  Copyright (C) 2000-2001 IDEALX
4%%%
5%%%  This program is free software; you can redistribute it and/or modify
6%%%  it under the terms of the GNU General Public License as published by
7%%%  the Free Software Foundation; either version 2 of the License, or
8%%%  (at your option) any later version.
9%%%
10%%%  This program is distributed in the hope that it will be useful,
11%%%  but WITHOUT ANY WARRANTY; without even the implied warranty of
12%%%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13%%%  GNU General Public License for more details.
14%%%
15%%%  You should have received a copy of the GNU General Public License
16%%%  along with this program; if not, write to the Free Software
17%%%  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
18%%%
19
20%%%  In addition, as a special exception, you have the permission to
21%%%  link the code of this program with any library released under
22%%%  the EPL license and distribute linked combinations including
23%%%  the two; the MPL (Mozilla Public License), which EPL (Erlang
24%%%  Public License) is based on, is included in this exception.
25
26%%% Random Generators for several probability distributions
27
28-module(ts_stats).
29-created('Date: 2000/10/20 13:58:56 nniclausse Exp ').
30-vc('$Id$ ').
31-author('nicolas.niclausse@niclux.org').
32
33-export([exponential/1, exponential/2, pareto/2,
34         normal/0, normal/1, normal/2, uniform/2,
35         invgaussian/2,
36         mean/1, mean/3,
37         variance/1,
38         meanvar/4,
39         meanvar_minmax/6,
40         stdvar/1]).
41
42-import(math, [log/1, pi/0, sqrt/1, pow/2]).
43
44-record(pareto, {a = 1 , beta}).
45-record(normal, {mean = 0 , stddev= 1 }).
46-record(invgaussian, {mu , lambda}).
47
48%% get n samples from a function F with parameter Param
49sample (F, Param, N) ->
50    sample(F, [], Param, N-1).
51
52sample (F, X, Param, 0) ->
53    [F(Param) | X] ;
54sample (F, X, Param, N) ->
55    sample(F, [F(Param)|X], Param, N-1 ).
56
57uniform(Min,Max)->
58    Min+random:uniform(Max-Min+1)-1.
59
60%% random sample from an exponential distribution
61exponential(Param) ->
62    -math:log(random:uniform())/Param.
63
64%% N samples from an exponential distribution
65exponential(Param, N) ->
66    sample(fun(X) -> exponential(X) end , Param, N).
67
68%% random sample from a Pareto distribution
69pareto(#pareto{a=A, beta=Beta}) ->
70    A/(math:pow(random:uniform(), 1/Beta)).
71
72%% if a list is given, construct a record for the parameters
73pareto([A, Beta], N) ->
74    pareto(#pareto{a = A , beta = Beta }, N);
75%% N samples from a Pareto distribution
76pareto(Param, N) ->
77    sample(fun(X) -> pareto(X) end , Param, N).
78
79invgaussian([Mu,Lambda],N) ->
80    invgaussian(#invgaussian{mu=Mu,lambda=Lambda},N);
81
82invgaussian(Param,N) ->
83    sample(fun(X) -> invgaussian(X) end , Param, N).
84
85%% random sample from a Inverse Gaussian distribution
86invgaussian(#invgaussian{mu=Mu, lambda=Lambda}) ->
87    Y = Mu*pow(normal(), 2),
88    X1 = Mu+Mu*Y/(2*Lambda)-Mu*sqrt(4*Lambda*Y+pow(Y,2))/(2*Lambda),
89    U = random:uniform(),
90    X = (Mu/(Mu+X1))-U,
91    case X >=0 of
92        true  -> X1;
93        false -> Mu*Mu/X1
94    end.
95
96normal() ->
97    [Val] = normal(#normal{},1),
98    Val.
99
100normal([Mean,StdDev],N) ->
101    normal(#normal{mean=Mean,stddev=StdDev},N);
102
103normal(Param,N) ->
104    sample(fun(X) -> normal(X) end , Param, N).
105
106normal(N) when is_integer(N)->
107    normal(#normal{},N);
108normal(#normal{mean=M,stddev=S}) ->
109    normal_boxm(M,S,0,0,1).
110
111%%% use the polar form of the Box-Muller transformation
112normal_boxm(M,S,X1,_X2,W) when W < 1->
113    W2 = sqrt( (-2.0 * log( W ) ) / W ),
114    Y1 = X1 * W2,
115    M + Y1 * S;
116normal_boxm(M,S,_,_,_W) ->
117    X1 = 2.0 * random:uniform() - 1.0,
118    X2 = 2.0 * random:uniform() - 1.0,
119    normal_boxm(M,S,X1,X2,X1 * X1 + X2 * X2).
120%%%
121
122%% incremental computation of the mean
123mean(Esp, [], _) -> Esp;
124
125mean(Esp, [X|H], I) ->
126    Next = I+1,
127    mean((Esp+(X-Esp)/(Next)), H, Next).
128
129%% compute the mean of a list
130mean([]) -> 0;
131
132mean(H) ->
133    mean(0, H, 0).
134
135%% @spec meanvar(Esp::number(),Var::number(),X::list() | number(),I::integer()) ->
136%%      {NewEsp::number(), NewVar::number(), Next::integer()}
137%% @doc incremental computation of the mean and variance together. The
138%%      algorithm should limit the round-off errors
139%% @end
140
141%% single value
142meanvar(Esp, Var, X, I) when is_number(X) ->
143    Next = I+1,
144    C = X - Esp,
145    NewEsp =  (X+Esp*I)/(Next),
146    NewVar = Var+C*(X-NewEsp),
147    { NewEsp, NewVar, Next };
148%% list of samples
149meanvar(Esp, Var,[], I) ->
150    {Esp, Var, I};
151meanvar(Esp, Var, [X|H], I) ->
152    {NewEsp, NewVar, Next} = meanvar(Esp,Var,X,I),
153    meanvar(NewEsp, NewVar, H, Next).
154
155%% compute min and max also
156meanvar_minmax(Esp, Var, Min, Max, X, I) when is_number(X)->
157    meanvar_minmax(Esp, Var, Min, Max, [X], I);
158meanvar_minmax(Esp, Var, Min, Max, [], I) ->
159    {Esp, Var, Min, Max, I};
160meanvar_minmax(Esp, Var, 0, 0, [X|H], I) -> % first data, set min and max
161    meanvar_minmax(Esp, Var, X, X, [X|H], I);
162meanvar_minmax(Esp, Var, Min, Max, [X|H], I) ->
163    {NewEsp, NewVar, Next} = meanvar(Esp,Var,X,I),
164    if
165        X > Max -> % new max, min unchanged
166            meanvar_minmax(NewEsp, NewVar, Min, X, H, Next);
167        X < Min -> % new min, max unchanged
168            meanvar_minmax(NewEsp, NewVar, X, Max, H, Next);
169        true ->
170            meanvar_minmax(NewEsp, NewVar, Min, Max, H, Next)
171    end.
172
173
174%% compute the variance of a list
175variance([]) -> 0;
176variance(H) ->
177    {_Mean, Var, I} = meanvar(0, 0, H, 0),
178    Var/I.
179
180stdvar(H) ->
181    math:sqrt(variance(H)).
182
183
184