1 // Copyright 2017 The Abseil Authors.
2 //
3 // Licensed under the Apache License, Version 2.0 (the "License");
4 // you may not use this file except in compliance with the License.
5 // You may obtain a copy of the License at
6 //
7 // https://www.apache.org/licenses/LICENSE-2.0
8 //
9 // Unless required by applicable law or agreed to in writing, software
10 // distributed under the License is distributed on an "AS IS" BASIS,
11 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12 // See the License for the specific language governing permissions and
13 // limitations under the License.
14
15 #include "absl/random/bernoulli_distribution.h"
16
17 #include <cmath>
18 #include <cstddef>
19 #include <random>
20 #include <sstream>
21 #include <utility>
22
23 #include "gtest/gtest.h"
24 #include "absl/random/internal/sequence_urbg.h"
25 #include "absl/random/random.h"
26
27 namespace {
28
29 class BernoulliTest : public testing::TestWithParam<std::pair<double, size_t>> {
30 };
31
TEST_P(BernoulliTest,Serialize)32 TEST_P(BernoulliTest, Serialize) {
33 const double d = GetParam().first;
34 absl::bernoulli_distribution before(d);
35
36 {
37 absl::bernoulli_distribution via_param{
38 absl::bernoulli_distribution::param_type(d)};
39 EXPECT_EQ(via_param, before);
40 }
41
42 std::stringstream ss;
43 ss << before;
44 absl::bernoulli_distribution after(0.6789);
45
46 EXPECT_NE(before.p(), after.p());
47 EXPECT_NE(before.param(), after.param());
48 EXPECT_NE(before, after);
49
50 ss >> after;
51
52 EXPECT_EQ(before.p(), after.p());
53 EXPECT_EQ(before.param(), after.param());
54 EXPECT_EQ(before, after);
55 }
56
TEST_P(BernoulliTest,Accuracy)57 TEST_P(BernoulliTest, Accuracy) {
58 // Sadly, the claim to fame for this implementation is precise accuracy, which
59 // is very, very hard to measure, the improvements come as trials approach the
60 // limit of double accuracy; thus the outcome differs from the
61 // std::bernoulli_distribution with a probability of approximately 1 in 2^-53.
62 const std::pair<double, size_t> para = GetParam();
63 size_t trials = para.second;
64 double p = para.first;
65
66 absl::InsecureBitGen rng;
67
68 size_t yes = 0;
69 absl::bernoulli_distribution dist(p);
70 for (size_t i = 0; i < trials; ++i) {
71 if (dist(rng)) yes++;
72 }
73
74 // Compute the distribution parameters for a binomial test, using a normal
75 // approximation for the confidence interval, as there are a sufficiently
76 // large number of trials that the central limit theorem applies.
77 const double stddev_p = std::sqrt((p * (1.0 - p)) / trials);
78 const double expected = trials * p;
79 const double stddev = trials * stddev_p;
80
81 // 5 sigma, approved by Richard Feynman
82 EXPECT_NEAR(yes, expected, 5 * stddev)
83 << "@" << p << ", "
84 << std::abs(static_cast<double>(yes) - expected) / stddev << " stddev";
85 }
86
87 // There must be many more trials to make the mean approximately normal for `p`
88 // closes to 0 or 1.
89 INSTANTIATE_TEST_SUITE_P(
90 All, BernoulliTest,
91 ::testing::Values(
92 // Typical values.
93 std::make_pair(0, 30000), std::make_pair(1e-3, 30000000),
94 std::make_pair(0.1, 3000000), std::make_pair(0.5, 3000000),
95 std::make_pair(0.9, 30000000), std::make_pair(0.999, 30000000),
96 std::make_pair(1, 30000),
97 // Boundary cases.
98 std::make_pair(std::nextafter(1.0, 0.0), 1), // ~1 - epsilon
99 std::make_pair(std::numeric_limits<double>::epsilon(), 1),
100 std::make_pair(std::nextafter(std::numeric_limits<double>::min(),
101 1.0), // min + epsilon
102 1),
103 std::make_pair(std::numeric_limits<double>::min(), // smallest normal
104 1),
105 std::make_pair(
106 std::numeric_limits<double>::denorm_min(), // smallest denorm
107 1),
108 std::make_pair(std::numeric_limits<double>::min() / 2, 1), // denorm
109 std::make_pair(std::nextafter(std::numeric_limits<double>::min(),
110 0.0), // denorm_max
111 1)));
112
113 // NOTE: absl::bernoulli_distribution is not guaranteed to be stable.
TEST(BernoulliTest,StabilityTest)114 TEST(BernoulliTest, StabilityTest) {
115 // absl::bernoulli_distribution stability relies on FastUniformBits and
116 // integer arithmetic.
117 absl::random_internal::sequence_urbg urbg({
118 0x0003eb76f6f7f755ull, 0xFFCEA50FDB2F953Bull, 0xC332DDEFBE6C5AA5ull,
119 0x6558218568AB9702ull, 0x2AEF7DAD5B6E2F84ull, 0x1521B62829076170ull,
120 0xECDD4775619F1510ull, 0x13CCA830EB61BD96ull, 0x0334FE1EAA0363CFull,
121 0xB5735C904C70A239ull, 0xD59E9E0BCBAADE14ull, 0xEECC86BC60622CA7ull,
122 0x4864f22c059bf29eull, 0x247856d8b862665cull, 0xe46e86e9a1337e10ull,
123 0xd8c8541f3519b133ull, 0xe75b5162c567b9e4ull, 0xf732e5ded7009c5bull,
124 0xb170b98353121eacull, 0x1ec2e8986d2362caull, 0x814c8e35fe9a961aull,
125 0x0c3cd59c9b638a02ull, 0xcb3bb6478a07715cull, 0x1224e62c978bbc7full,
126 0x671ef2cb04e81f6eull, 0x3c1cbd811eaf1808ull, 0x1bbc23cfa8fac721ull,
127 0xa4c2cda65e596a51ull, 0xb77216fad37adf91ull, 0x836d794457c08849ull,
128 0xe083df03475f49d7ull, 0xbc9feb512e6b0d6cull, 0xb12d74fdd718c8c5ull,
129 0x12ff09653bfbe4caull, 0x8dd03a105bc4ee7eull, 0x5738341045ba0d85ull,
130 0xe3fd722dc65ad09eull, 0x5a14fd21ea2a5705ull, 0x14e6ea4d6edb0c73ull,
131 0x275b0dc7e0a18acfull, 0x36cebe0d2653682eull, 0x0361e9b23861596bull,
132 });
133
134 // Generate a std::string of '0' and '1' for the distribution output.
135 auto generate = [&urbg](absl::bernoulli_distribution& dist) {
136 std::string output;
137 output.reserve(36);
138 urbg.reset();
139 for (int i = 0; i < 35; i++) {
140 output.append(dist(urbg) ? "1" : "0");
141 }
142 return output;
143 };
144
145 const double kP = 0.0331289862362;
146 {
147 absl::bernoulli_distribution dist(kP);
148 auto v = generate(dist);
149 EXPECT_EQ(35, urbg.invocations());
150 EXPECT_EQ(v, "00000000000010000000000010000000000") << dist;
151 }
152 {
153 absl::bernoulli_distribution dist(kP * 10.0);
154 auto v = generate(dist);
155 EXPECT_EQ(35, urbg.invocations());
156 EXPECT_EQ(v, "00000100010010010010000011000011010") << dist;
157 }
158 {
159 absl::bernoulli_distribution dist(kP * 20.0);
160 auto v = generate(dist);
161 EXPECT_EQ(35, urbg.invocations());
162 EXPECT_EQ(v, "00011110010110110011011111110111011") << dist;
163 }
164 {
165 absl::bernoulli_distribution dist(1.0 - kP);
166 auto v = generate(dist);
167 EXPECT_EQ(35, urbg.invocations());
168 EXPECT_EQ(v, "11111111111111111111011111111111111") << dist;
169 }
170 }
171
TEST(BernoulliTest,StabilityTest2)172 TEST(BernoulliTest, StabilityTest2) {
173 absl::random_internal::sequence_urbg urbg(
174 {0x0003eb76f6f7f755ull, 0xFFCEA50FDB2F953Bull, 0xC332DDEFBE6C5AA5ull,
175 0x6558218568AB9702ull, 0x2AEF7DAD5B6E2F84ull, 0x1521B62829076170ull,
176 0xECDD4775619F1510ull, 0x13CCA830EB61BD96ull, 0x0334FE1EAA0363CFull,
177 0xB5735C904C70A239ull, 0xD59E9E0BCBAADE14ull, 0xEECC86BC60622CA7ull});
178
179 // Generate a std::string of '0' and '1' for the distribution output.
180 auto generate = [&urbg](absl::bernoulli_distribution& dist) {
181 std::string output;
182 output.reserve(13);
183 urbg.reset();
184 for (int i = 0; i < 12; i++) {
185 output.append(dist(urbg) ? "1" : "0");
186 }
187 return output;
188 };
189
190 constexpr double b0 = 1.0 / 13.0 / 0.2;
191 constexpr double b1 = 2.0 / 13.0 / 0.2;
192 constexpr double b3 = (5.0 / 13.0 / 0.2) - ((1 - b0) + (1 - b1) + (1 - b1));
193 {
194 absl::bernoulli_distribution dist(b0);
195 auto v = generate(dist);
196 EXPECT_EQ(12, urbg.invocations());
197 EXPECT_EQ(v, "000011100101") << dist;
198 }
199 {
200 absl::bernoulli_distribution dist(b1);
201 auto v = generate(dist);
202 EXPECT_EQ(12, urbg.invocations());
203 EXPECT_EQ(v, "001111101101") << dist;
204 }
205 {
206 absl::bernoulli_distribution dist(b3);
207 auto v = generate(dist);
208 EXPECT_EQ(12, urbg.invocations());
209 EXPECT_EQ(v, "001111101111") << dist;
210 }
211 }
212
213 } // namespace
214