Adding random variable generation functions

src/stats_rv.erl

@@ -0,0 +1,73 @@
+%% -------------------------------------------------------------------
+%%
+%% stats: Statistics Suite for Erlang
+%%
+%% Copyright (c) 2009 Dave Smith (dizzyd@dizzyd.com)
+%%
+%% This file is provided to you under the Apache License,
+%% Version 2.0 (the "License"); you may not use this file
+%% except in compliance with the License.  You may obtain
+%% a copy of the License at
+%%
+%%   http://www.apache.org/licenses/LICENSE-2.0
+%%
+%% Unless required by applicable law or agreed to in writing,
+%% software distributed under the License is distributed on an
+%% "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+%% KIND, either express or implied.  See the License for the
+%% specific language governing permissions and limitations
+%% under the License.
+%%
+%% -------------------------------------------------------------------
+-module(stats_rv).
+
+-export([uniform/0,
+         exponential/1,
+         poisson/1,
+         normal/2]).
+
+
+%% ====================================================================
+%% Public API
+%% ====================================================================
+
+%%
+%% Generates a uniformly-distributed random variable (wrapper for convenience)
+%%
+uniform() ->
+    random:uniform().
+
+%%
+%% Generates an exponential-distributed random variable, using inverse function
+%%
+exponential(Lambda) ->
+    -math:log(random:uniform()) / Lambda.
+
+%%
+%% Generates a Poisson-distributed random variable, using Knuth's method
+%%
+poisson(Lambda) ->
+    poisson_rv_loop(math:exp(-Lambda), 1, 0).
+
+%%
+%% Generates a Normal-distributed random variable, using Box-Muller method
+%%
+normal(Mean, Sigma) ->
+    Rv1 = random:uniform(),
+    Rv2 = random:uniform(),
+    Rho = math:sqrt(-2 * math:log(1-Rv2)),
+    Rho * math:cos(2 * math:pi() * Rv1) * Sigma + Mean.
+
+
+
+%% ====================================================================
+%% Internal functions
+%% ====================================================================
+
+poisson_rv_loop(A, B, Counter) ->
+    case B < A of
+        true ->
+            Counter;
+        false ->
+            poisson_rv_loop(A, B * random:uniform(), Counter+1)
+    end.