Test normal distribution

1 files changed, 204 insertions, 1 deletions

diff --git a/lib/stdlib/test/rand_SUITE.erl b/lib/stdlib/test/rand_SUITE.erl
index ef4f9faad9..dbf1857a25 100644
--- a/lib/stdlib/test/rand_SUITE.erl
+++ b/lib/stdlib/test/rand_SUITE.erl
@@ -29,6 +29,9 @@
 	 basic_stats_uniform_1/1, basic_stats_uniform_2/1,
 	 basic_stats_standard_normal/1,
 	 basic_stats_normal/1,
+         stats_standard_normal_box_muller/1,
+         stats_standard_normal_box_muller_2/1,
+         stats_standard_normal/1,
          uniform_real_conv/1,
 	 plugin/1, measure/1,
 	 reference_jump_state/1, reference_jump_procdict/1]).
@@ -57,7 +60,10 @@ all() ->
 groups() ->
     [{basic_stats, [parallel],
       [basic_stats_uniform_1, basic_stats_uniform_2,
-       basic_stats_standard_normal]},
+       basic_stats_standard_normal,
+       stats_standard_normal_box_muller,
+       stats_standard_normal_box_muller_2,
+       stats_standard_normal]},
      {reference_jump, [parallel],
       [reference_jump_state, reference_jump_procdict]}].
 
@@ -410,6 +416,203 @@ normal_s(Mean, Variance, State0) when Mean == 0, Variance == 1 ->
 normal_s(Mean, Variance, State0) ->
     rand:normal_s(Mean, Variance, State0).
 
+
+
+-dialyzer({no_improper_lists, stats_standard_normal_box_muller/1}).
+stats_standard_normal_box_muller(Config) when is_list(Config) ->
+    try math:erfc(1.0) of
+        _ ->
+            TwoPi = 2.0 * math:pi(),
+            NormalS =
+                fun
+                    ([S0]) ->
+                        {U1, S1} = rand:uniform_real_s(S0),
+                        R = math:sqrt(-2.0 * math:log(U1)),
+                        {U2, S2} = rand:uniform_s(S1),
+                        T = TwoPi * U2,
+                        Z0 = R * math:cos(T),
+                        Z1 = R * math:sin(T),
+                        {Z0, [S2|Z1]};
+                    ([S|Z]) ->
+                        {Z, [S]}
+                end,
+            State = [rand:seed(exrop)],
+            stats_standard_normal(NormalS, State)
+    catch error:_ ->
+            {skip, "math:erfc/1 not supported"}
+    end.
+
+-dialyzer({no_improper_lists, stats_standard_normal_box_muller_2/1}).
+stats_standard_normal_box_muller_2(Config) when is_list(Config) ->
+    try math:erfc(1.0) of
+        _ ->
+            TwoPi = 2.0 * math:pi(),
+            NormalS =
+                fun
+                    ([S0]) ->
+                        {U0, S1} = rand:uniform_s(S0),
+                        U1 = 1.0 - U0,
+                        R = math:sqrt(-2.0 * math:log(U1)),
+                        {U2, S2} = rand:uniform_s(S1),
+                        T = TwoPi * U2,
+                        Z0 = R * math:cos(T),
+                        Z1 = R * math:sin(T),
+                        {Z0, [S2|Z1]};
+                    ([S|Z]) ->
+                        {Z, [S]}
+                end,
+            State = [rand:seed(exrop)],
+            stats_standard_normal(NormalS, State)
+    catch error:_ ->
+            {skip, "math:erfc/1 not supported"}
+    end.
+
+
+stats_standard_normal(Config) when is_list(Config) ->
+    try math:erfc(1.0) of
+        _ ->
+            stats_standard_normal(
+              fun rand:normal_s/1, rand:seed_s(exrop))
+    catch error:_ ->
+            {skip, "math:erfc/1 not supported"}
+    end.
+%%
+stats_standard_normal(Fun, S) ->
+%%%
+%%% ct config:
+%%% {rand_SUITE, [{stats_standard_normal,[{seconds, 8}, {std_devs, 4.0}]}]}.
+%%%
+    Seconds = ct:get_config({?MODULE, ?FUNCTION_NAME, seconds}, 8),
+    StdDevs =
+        ct:get_config(
+          {?MODULE, ?FUNCTION_NAME, std_devs},
+          4.0), % probability erfc(4/sqrt(2)) (1/15787) to fail a bucket
+%%%
+    ct:timetrap({seconds, Seconds + 120}),
+    %% InvDelta and Buckets are chosen so that the probability to land
+    %% in the central bucket (closest to 0) is about the same as
+    %% the probability to land in the top catch-all-outside bucket.
+    %%
+    %% Rounds is calculated so the expected value for the central
+    %% bucket will be at least TargetHits.
+    %%
+    InvDelta = 512,
+    W = InvDelta * math:sqrt(2),
+    Buckets = 1619, % erf-1(erfc(1 / W)) * W
+    TargetHits = 1024,
+    InvP0 = floor(math:ceil(1 / (0.5 * math:erf(1 / W)))),
+    Rounds = 2 * TargetHits * InvP0,
+    Histogram = array:new({default, 0}),
+    StopTime = erlang:monotonic_time(second) + Seconds,
+    ct:pal(
+      "Running standard normal test against ~w std devs for ~w seconds...",
+      [StdDevs, Seconds]),
+    {PositiveHistogram, NegativeHistogram, Outlier, TotalRounds} =
+        stats_standard_normal(
+          InvDelta, Buckets, Histogram, Histogram, 0.0,
+          Fun, S, Rounds, StopTime, Rounds, 0),
+    Precision = math:sqrt(TotalRounds / InvP0) / StdDevs,
+    ct:pal(
+      "Total rounds: ~w, tolerance: 1/~.2f, outlier: ~.2f.",
+      [TotalRounds, Precision, Outlier]),
+    ResultingTargetHits = floor(TotalRounds / 2 / InvP0),
+    {TotalRounds, [], []} =
+        {TotalRounds,
+         check_histogram(
+           W, Buckets, ResultingTargetHits, TotalRounds,
+           StdDevs, PositiveHistogram),
+         check_histogram(
+           W, Buckets, ResultingTargetHits, TotalRounds,
+           StdDevs, NegativeHistogram)},
+    ok.
+%%
+stats_standard_normal(
+  InvDelta, Buckets, PositiveHistogram, NegativeHistogram, Outlier,
+  Fun, S, 0, StopTime, Rounds, TotalRounds) ->
+    case erlang:monotonic_time(second) of
+        Now when Now < StopTime ->
+            stats_standard_normal(
+              InvDelta, Buckets,
+              PositiveHistogram, NegativeHistogram, Outlier,
+              Fun, S, Rounds, StopTime, Rounds, TotalRounds + Rounds);
+        _ ->
+            {PositiveHistogram, NegativeHistogram,
+             Outlier, TotalRounds + Rounds}
+    end;
+stats_standard_normal(
+  InvDelta, Buckets, PositiveHistogram, NegativeHistogram, Outlier,
+  Fun, S, Count, StopTime, Rounds, TotalRounds) ->
+    case Fun(S) of
+        {X, NewS} when 0.0 =< X ->
+            Bucket = min(Buckets, floor(X * InvDelta)),
+            stats_standard_normal(
+              InvDelta, Buckets,
+              increment_bucket(Bucket, PositiveHistogram),
+              NegativeHistogram, max(Outlier, X),
+              Fun, NewS, Count - 1, StopTime, Rounds, TotalRounds);
+        {MinusX, NewS} ->
+            X = -MinusX,
+            Bucket = min(Buckets, floor(X * InvDelta)),
+            stats_standard_normal(
+              InvDelta, Buckets,
+              PositiveHistogram,
+              increment_bucket(Bucket, NegativeHistogram), max(Outlier, X),
+              Fun, NewS, Count - 1, StopTime, Rounds, TotalRounds)
+    end.
+
+increment_bucket(Bucket, Array) ->
+    array:set(Bucket, array:get(Bucket, Array) + 1, Array).
+
+check_histogram(W, Buckets, TargetHits, Rounds, StdDevs, Histogram) ->
+    %%PrevBucket = 512,
+    %%Bucket = PrevBucket - 1,
+    %%P = 0.5 * math:erfc(PrevBucket / W),
+    PrevBucket = Buckets,
+    Bucket = PrevBucket,
+    P = 0.0,
+    N = 0,
+    check_histogram(
+      W, Bucket, TargetHits, Rounds, StdDevs, Histogram,
+      PrevBucket, P, N).
+%%
+check_histogram(
+  _W, _Bucket, _TargetHits, _Rounds, _StdDevs, _Histogram,
+  0, _PrevP, _PrevN) ->
+    [];
+check_histogram(
+  W, Bucket, TargetHits, Rounds, StdDevs, Histogram,
+  PrevBucket, PrevP, PrevN) ->
+    N = PrevN + array:get(Bucket, Histogram),
+    P = 0.5 * math:erfc(Bucket / W),
+    DP = P - PrevP,
+    Exp = DP * Rounds,
+    if
+        TargetHits =< Exp ->
+            %% Check N
+            Err = N - Exp,
+            AbsErr = abs(Err),
+            Var = Rounds * DP*(1.0 - DP),
+            Limit = StdDevs * math:sqrt(Var),
+            if
+                Limit =< AbsErr ->
+                    [#{bucket => {Bucket, PrevBucket},
+                       n => N, exp => Exp, err => Err, limit => Limit} |
+                     check_histogram(
+                       W, Bucket - 1, TargetHits, Rounds, StdDevs, Histogram,
+                       Bucket, P, 0)];
+                    %% ct:fail(
+                    %%   {bucket, Bucket, PrevBucket, N, Exp, Err, Limit});
+                true ->
+                    check_histogram(
+                      W, Bucket - 1, TargetHits, Rounds, StdDevs, Histogram,
+                      Bucket, P, 0)
+            end;
+        true ->
+            check_histogram(
+              W, Bucket - 1, TargetHits, Rounds, StdDevs, Histogram,
+              PrevBucket, PrevP, N)
+    end.
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% White box test of the conversion to float