aboutsummaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorRaimo Niskanen <[email protected]>2017-10-22 21:01:06 +0200
committerRaimo Niskanen <[email protected]>2017-10-27 11:54:22 +0200
commit15a942b0f2c68ce8a0398db37adb020fd4063b58 (patch)
tree669341ba54b4e9013bec4f00c18f1a3499ca2306
parent5ce0138c0809bd3f17029413fdf2ead1a8979762 (diff)
downloadotp-15a942b0f2c68ce8a0398db37adb020fd4063b58.tar.gz
otp-15a942b0f2c68ce8a0398db37adb020fd4063b58.tar.bz2
otp-15a942b0f2c68ce8a0398db37adb020fd4063b58.zip
Test normal distribution
-rw-r--r--lib/stdlib/test/rand_SUITE.erl205
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