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