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/beta_distribution.h"
16
17 #include <algorithm>
18 #include <cfloat>
19 #include <cstddef>
20 #include <cstdint>
21 #include <iterator>
22 #include <random>
23 #include <sstream>
24 #include <string>
25 #include <type_traits>
26 #include <unordered_map>
27 #include <vector>
28
29 #include "gmock/gmock.h"
30 #include "gtest/gtest.h"
31 #include "absl/base/internal/raw_logging.h"
32 #include "absl/numeric/internal/representation.h"
33 #include "absl/random/internal/chi_square.h"
34 #include "absl/random/internal/distribution_test_util.h"
35 #include "absl/random/internal/pcg_engine.h"
36 #include "absl/random/internal/sequence_urbg.h"
37 #include "absl/random/random.h"
38 #include "absl/strings/str_cat.h"
39 #include "absl/strings/str_format.h"
40 #include "absl/strings/str_replace.h"
41 #include "absl/strings/strip.h"
42
43 namespace {
44
45 template <typename IntType>
46 class BetaDistributionInterfaceTest : public ::testing::Test {};
47
ShouldExerciseLongDoubleTests()48 constexpr bool ShouldExerciseLongDoubleTests() {
49 // long double arithmetic is not supported well by either GCC or Clang on
50 // most platforms specifically not when implemented in terms of double-double;
51 // see https://gcc.gnu.org/bugzilla/show_bug.cgi?id=99048,
52 // https://bugs.llvm.org/show_bug.cgi?id=49131, and
53 // https://bugs.llvm.org/show_bug.cgi?id=49132.
54 // So a conservative choice here is to disable long-double tests pretty much
55 // everywhere except on x64 but only if long double is not implemented as
56 // double-double.
57 #if defined(__i686__) && defined(__x86_64__)
58 return !absl::numeric_internal::IsDoubleDouble();
59 #else
60 return false;
61 #endif
62 }
63
64 using RealTypes = std::conditional<ShouldExerciseLongDoubleTests(),
65 ::testing::Types<float, double, long double>,
66 ::testing::Types<float, double>>::type;
67 TYPED_TEST_SUITE(BetaDistributionInterfaceTest, RealTypes);
68
TYPED_TEST(BetaDistributionInterfaceTest,SerializeTest)69 TYPED_TEST(BetaDistributionInterfaceTest, SerializeTest) {
70 // The threshold for whether std::exp(1/a) is finite.
71 const TypeParam kSmallA =
72 1.0f / std::log((std::numeric_limits<TypeParam>::max)());
73 // The threshold for whether a * std::log(a) is finite.
74 const TypeParam kLargeA =
75 std::exp(std::log((std::numeric_limits<TypeParam>::max)()) -
76 std::log(std::log((std::numeric_limits<TypeParam>::max)())));
77 using param_type = typename absl::beta_distribution<TypeParam>::param_type;
78
79 constexpr int kCount = 1000;
80 absl::InsecureBitGen gen;
81 const TypeParam kValues[] = {
82 TypeParam(1e-20), TypeParam(1e-12), TypeParam(1e-8), TypeParam(1e-4),
83 TypeParam(1e-3), TypeParam(0.1), TypeParam(0.25),
84 std::nextafter(TypeParam(0.5), TypeParam(0)), // 0.5 - epsilon
85 std::nextafter(TypeParam(0.5), TypeParam(1)), // 0.5 + epsilon
86 TypeParam(0.5), TypeParam(1.0), //
87 std::nextafter(TypeParam(1), TypeParam(0)), // 1 - epsilon
88 std::nextafter(TypeParam(1), TypeParam(2)), // 1 + epsilon
89 TypeParam(12.5), TypeParam(1e2), TypeParam(1e8), TypeParam(1e12),
90 TypeParam(1e20), //
91 kSmallA, //
92 std::nextafter(kSmallA, TypeParam(0)), //
93 std::nextafter(kSmallA, TypeParam(1)), //
94 kLargeA, //
95 std::nextafter(kLargeA, TypeParam(0)), //
96 std::nextafter(kLargeA, std::numeric_limits<TypeParam>::max()),
97 // Boundary cases.
98 std::numeric_limits<TypeParam>::max(),
99 std::numeric_limits<TypeParam>::epsilon(),
100 std::nextafter(std::numeric_limits<TypeParam>::min(),
101 TypeParam(1)), // min + epsilon
102 std::numeric_limits<TypeParam>::min(), // smallest normal
103 std::numeric_limits<TypeParam>::denorm_min(), // smallest denorm
104 std::numeric_limits<TypeParam>::min() / 2, // denorm
105 std::nextafter(std::numeric_limits<TypeParam>::min(),
106 TypeParam(0)), // denorm_max
107 };
108 for (TypeParam alpha : kValues) {
109 for (TypeParam beta : kValues) {
110 ABSL_INTERNAL_LOG(
111 INFO, absl::StrFormat("Smoke test for Beta(%a, %a)", alpha, beta));
112
113 param_type param(alpha, beta);
114 absl::beta_distribution<TypeParam> before(alpha, beta);
115 EXPECT_EQ(before.alpha(), param.alpha());
116 EXPECT_EQ(before.beta(), param.beta());
117
118 {
119 absl::beta_distribution<TypeParam> via_param(param);
120 EXPECT_EQ(via_param, before);
121 EXPECT_EQ(via_param.param(), before.param());
122 }
123
124 // Smoke test.
125 for (int i = 0; i < kCount; ++i) {
126 auto sample = before(gen);
127 EXPECT_TRUE(std::isfinite(sample));
128 EXPECT_GE(sample, before.min());
129 EXPECT_LE(sample, before.max());
130 }
131
132 // Validate stream serialization.
133 std::stringstream ss;
134 ss << before;
135 absl::beta_distribution<TypeParam> after(3.8f, 1.43f);
136 EXPECT_NE(before.alpha(), after.alpha());
137 EXPECT_NE(before.beta(), after.beta());
138 EXPECT_NE(before.param(), after.param());
139 EXPECT_NE(before, after);
140
141 ss >> after;
142
143 EXPECT_EQ(before.alpha(), after.alpha());
144 EXPECT_EQ(before.beta(), after.beta());
145 EXPECT_EQ(before, after) //
146 << ss.str() << " " //
147 << (ss.good() ? "good " : "") //
148 << (ss.bad() ? "bad " : "") //
149 << (ss.eof() ? "eof " : "") //
150 << (ss.fail() ? "fail " : "");
151 }
152 }
153 }
154
TYPED_TEST(BetaDistributionInterfaceTest,DegenerateCases)155 TYPED_TEST(BetaDistributionInterfaceTest, DegenerateCases) {
156 // We use a fixed bit generator for distribution accuracy tests. This allows
157 // these tests to be deterministic, while still testing the qualify of the
158 // implementation.
159 absl::random_internal::pcg64_2018_engine rng(0x2B7E151628AED2A6);
160
161 // Extreme cases when the params are abnormal.
162 constexpr int kCount = 1000;
163 const TypeParam kSmallValues[] = {
164 std::numeric_limits<TypeParam>::min(),
165 std::numeric_limits<TypeParam>::denorm_min(),
166 std::nextafter(std::numeric_limits<TypeParam>::min(),
167 TypeParam(0)), // denorm_max
168 std::numeric_limits<TypeParam>::epsilon(),
169 };
170 const TypeParam kLargeValues[] = {
171 std::numeric_limits<TypeParam>::max() * static_cast<TypeParam>(0.9999),
172 std::numeric_limits<TypeParam>::max() - 1,
173 std::numeric_limits<TypeParam>::max(),
174 };
175 {
176 // Small alpha and beta.
177 // Useful WolframAlpha plots:
178 // * plot InverseBetaRegularized[x, 0.0001, 0.0001] from 0.495 to 0.505
179 // * Beta[1.0, 0.0000001, 0.0000001]
180 // * Beta[0.9999, 0.0000001, 0.0000001]
181 for (TypeParam alpha : kSmallValues) {
182 for (TypeParam beta : kSmallValues) {
183 int zeros = 0;
184 int ones = 0;
185 absl::beta_distribution<TypeParam> d(alpha, beta);
186 for (int i = 0; i < kCount; ++i) {
187 TypeParam x = d(rng);
188 if (x == 0.0) {
189 zeros++;
190 } else if (x == 1.0) {
191 ones++;
192 }
193 }
194 EXPECT_EQ(ones + zeros, kCount);
195 if (alpha == beta) {
196 EXPECT_NE(ones, 0);
197 EXPECT_NE(zeros, 0);
198 }
199 }
200 }
201 }
202 {
203 // Small alpha, large beta.
204 // Useful WolframAlpha plots:
205 // * plot InverseBetaRegularized[x, 0.0001, 10000] from 0.995 to 1
206 // * Beta[0, 0.0000001, 1000000]
207 // * Beta[0.001, 0.0000001, 1000000]
208 // * Beta[1, 0.0000001, 1000000]
209 for (TypeParam alpha : kSmallValues) {
210 for (TypeParam beta : kLargeValues) {
211 absl::beta_distribution<TypeParam> d(alpha, beta);
212 for (int i = 0; i < kCount; ++i) {
213 EXPECT_EQ(d(rng), 0.0);
214 }
215 }
216 }
217 }
218 {
219 // Large alpha, small beta.
220 // Useful WolframAlpha plots:
221 // * plot InverseBetaRegularized[x, 10000, 0.0001] from 0 to 0.001
222 // * Beta[0.99, 1000000, 0.0000001]
223 // * Beta[1, 1000000, 0.0000001]
224 for (TypeParam alpha : kLargeValues) {
225 for (TypeParam beta : kSmallValues) {
226 absl::beta_distribution<TypeParam> d(alpha, beta);
227 for (int i = 0; i < kCount; ++i) {
228 EXPECT_EQ(d(rng), 1.0);
229 }
230 }
231 }
232 }
233 {
234 // Large alpha and beta.
235 absl::beta_distribution<TypeParam> d(std::numeric_limits<TypeParam>::max(),
236 std::numeric_limits<TypeParam>::max());
237 for (int i = 0; i < kCount; ++i) {
238 EXPECT_EQ(d(rng), 0.5);
239 }
240 }
241 {
242 // Large alpha and beta but unequal.
243 absl::beta_distribution<TypeParam> d(
244 std::numeric_limits<TypeParam>::max(),
245 std::numeric_limits<TypeParam>::max() * 0.9999);
246 for (int i = 0; i < kCount; ++i) {
247 TypeParam x = d(rng);
248 EXPECT_NE(x, 0.5f);
249 EXPECT_FLOAT_EQ(x, 0.500025f);
250 }
251 }
252 }
253
254 class BetaDistributionModel {
255 public:
BetaDistributionModel(::testing::tuple<double,double> p)256 explicit BetaDistributionModel(::testing::tuple<double, double> p)
257 : alpha_(::testing::get<0>(p)), beta_(::testing::get<1>(p)) {}
258
Mean() const259 double Mean() const { return alpha_ / (alpha_ + beta_); }
260
Variance() const261 double Variance() const {
262 return alpha_ * beta_ / (alpha_ + beta_ + 1) / (alpha_ + beta_) /
263 (alpha_ + beta_);
264 }
265
Kurtosis() const266 double Kurtosis() const {
267 return 3 + 6 *
268 ((alpha_ - beta_) * (alpha_ - beta_) * (alpha_ + beta_ + 1) -
269 alpha_ * beta_ * (2 + alpha_ + beta_)) /
270 alpha_ / beta_ / (alpha_ + beta_ + 2) / (alpha_ + beta_ + 3);
271 }
272
273 protected:
274 const double alpha_;
275 const double beta_;
276 };
277
278 class BetaDistributionTest
279 : public ::testing::TestWithParam<::testing::tuple<double, double>>,
280 public BetaDistributionModel {
281 public:
BetaDistributionTest()282 BetaDistributionTest() : BetaDistributionModel(GetParam()) {}
283
284 protected:
285 template <class D>
286 bool SingleZTestOnMeanAndVariance(double p, size_t samples);
287
288 template <class D>
289 bool SingleChiSquaredTest(double p, size_t samples, size_t buckets);
290
291 absl::InsecureBitGen rng_;
292 };
293
294 template <class D>
SingleZTestOnMeanAndVariance(double p,size_t samples)295 bool BetaDistributionTest::SingleZTestOnMeanAndVariance(double p,
296 size_t samples) {
297 D dis(alpha_, beta_);
298
299 std::vector<double> data;
300 data.reserve(samples);
301 for (size_t i = 0; i < samples; i++) {
302 const double variate = dis(rng_);
303 EXPECT_FALSE(std::isnan(variate));
304 // Note that equality is allowed on both sides.
305 EXPECT_GE(variate, 0.0);
306 EXPECT_LE(variate, 1.0);
307 data.push_back(variate);
308 }
309
310 // We validate that the sample mean and sample variance are indeed from a
311 // Beta distribution with the given shape parameters.
312 const auto m = absl::random_internal::ComputeDistributionMoments(data);
313
314 // The variance of the sample mean is variance / n.
315 const double mean_stddev = std::sqrt(Variance() / static_cast<double>(m.n));
316
317 // The variance of the sample variance is (approximately):
318 // (kurtosis - 1) * variance^2 / n
319 const double variance_stddev = std::sqrt(
320 (Kurtosis() - 1) * Variance() * Variance() / static_cast<double>(m.n));
321 // z score for the sample variance.
322 const double z_variance = (m.variance - Variance()) / variance_stddev;
323
324 const double max_err = absl::random_internal::MaxErrorTolerance(p);
325 const double z_mean = absl::random_internal::ZScore(Mean(), m);
326 const bool pass =
327 absl::random_internal::Near("z", z_mean, 0.0, max_err) &&
328 absl::random_internal::Near("z_variance", z_variance, 0.0, max_err);
329 if (!pass) {
330 ABSL_INTERNAL_LOG(
331 INFO,
332 absl::StrFormat(
333 "Beta(%f, %f), "
334 "mean: sample %f, expect %f, which is %f stddevs away, "
335 "variance: sample %f, expect %f, which is %f stddevs away.",
336 alpha_, beta_, m.mean, Mean(),
337 std::abs(m.mean - Mean()) / mean_stddev, m.variance, Variance(),
338 std::abs(m.variance - Variance()) / variance_stddev));
339 }
340 return pass;
341 }
342
343 template <class D>
SingleChiSquaredTest(double p,size_t samples,size_t buckets)344 bool BetaDistributionTest::SingleChiSquaredTest(double p, size_t samples,
345 size_t buckets) {
346 constexpr double kErr = 1e-7;
347 std::vector<double> cutoffs, expected;
348 const double bucket_width = 1.0 / static_cast<double>(buckets);
349 int i = 1;
350 int unmerged_buckets = 0;
351 for (; i < buckets; ++i) {
352 const double p = bucket_width * static_cast<double>(i);
353 const double boundary =
354 absl::random_internal::BetaIncompleteInv(alpha_, beta_, p);
355 // The intention is to add `boundary` to the list of `cutoffs`. It becomes
356 // problematic, however, when the boundary values are not monotone, due to
357 // numerical issues when computing the inverse regularized incomplete
358 // Beta function. In these cases, we merge that bucket with its previous
359 // neighbor and merge their expected counts.
360 if ((cutoffs.empty() && boundary < kErr) ||
361 (!cutoffs.empty() && boundary <= cutoffs.back())) {
362 unmerged_buckets++;
363 continue;
364 }
365 if (boundary >= 1.0 - 1e-10) {
366 break;
367 }
368 cutoffs.push_back(boundary);
369 expected.push_back(static_cast<double>(1 + unmerged_buckets) *
370 bucket_width * static_cast<double>(samples));
371 unmerged_buckets = 0;
372 }
373 cutoffs.push_back(std::numeric_limits<double>::infinity());
374 // Merge all remaining buckets.
375 expected.push_back(static_cast<double>(buckets - i + 1) * bucket_width *
376 static_cast<double>(samples));
377 // Make sure that we don't merge all the buckets, making this test
378 // meaningless.
379 EXPECT_GE(cutoffs.size(), 3) << alpha_ << ", " << beta_;
380
381 D dis(alpha_, beta_);
382
383 std::vector<int32_t> counts(cutoffs.size(), 0);
384 for (int i = 0; i < samples; i++) {
385 const double x = dis(rng_);
386 auto it = std::upper_bound(cutoffs.begin(), cutoffs.end(), x);
387 counts[std::distance(cutoffs.begin(), it)]++;
388 }
389
390 // Null-hypothesis is that the distribution is beta distributed with the
391 // provided alpha, beta params (not estimated from the data).
392 const int dof = cutoffs.size() - 1;
393
394 const double chi_square = absl::random_internal::ChiSquare(
395 counts.begin(), counts.end(), expected.begin(), expected.end());
396 const bool pass =
397 (absl::random_internal::ChiSquarePValue(chi_square, dof) >= p);
398 if (!pass) {
399 for (int i = 0; i < cutoffs.size(); i++) {
400 ABSL_INTERNAL_LOG(
401 INFO, absl::StrFormat("cutoff[%d] = %f, actual count %d, expected %d",
402 i, cutoffs[i], counts[i],
403 static_cast<int>(expected[i])));
404 }
405
406 ABSL_INTERNAL_LOG(
407 INFO, absl::StrFormat(
408 "Beta(%f, %f) %s %f, p = %f", alpha_, beta_,
409 absl::random_internal::kChiSquared, chi_square,
410 absl::random_internal::ChiSquarePValue(chi_square, dof)));
411 }
412 return pass;
413 }
414
TEST_P(BetaDistributionTest,TestSampleStatistics)415 TEST_P(BetaDistributionTest, TestSampleStatistics) {
416 static constexpr int kRuns = 20;
417 static constexpr double kPFail = 0.02;
418 const double p =
419 absl::random_internal::RequiredSuccessProbability(kPFail, kRuns);
420 static constexpr int kSampleCount = 10000;
421 static constexpr int kBucketCount = 100;
422 int failed = 0;
423 for (int i = 0; i < kRuns; ++i) {
424 if (!SingleZTestOnMeanAndVariance<absl::beta_distribution<double>>(
425 p, kSampleCount)) {
426 failed++;
427 }
428 if (!SingleChiSquaredTest<absl::beta_distribution<double>>(
429 0.005, kSampleCount, kBucketCount)) {
430 failed++;
431 }
432 }
433 // Set so that the test is not flaky at --runs_per_test=10000
434 EXPECT_LE(failed, 5);
435 }
436
ParamName(const::testing::TestParamInfo<::testing::tuple<double,double>> & info)437 std::string ParamName(
438 const ::testing::TestParamInfo<::testing::tuple<double, double>>& info) {
439 std::string name = absl::StrCat("alpha_", ::testing::get<0>(info.param),
440 "__beta_", ::testing::get<1>(info.param));
441 return absl::StrReplaceAll(name, {{"+", "_"}, {"-", "_"}, {".", "_"}});
442 }
443
444 INSTANTIATE_TEST_SUITE_P(
445 TestSampleStatisticsCombinations, BetaDistributionTest,
446 ::testing::Combine(::testing::Values(0.1, 0.2, 0.9, 1.1, 2.5, 10.0, 123.4),
447 ::testing::Values(0.1, 0.2, 0.9, 1.1, 2.5, 10.0, 123.4)),
448 ParamName);
449
450 INSTANTIATE_TEST_SUITE_P(
451 TestSampleStatistics_SelectedPairs, BetaDistributionTest,
452 ::testing::Values(std::make_pair(0.5, 1000), std::make_pair(1000, 0.5),
453 std::make_pair(900, 1000), std::make_pair(10000, 20000),
454 std::make_pair(4e5, 2e7), std::make_pair(1e7, 1e5)),
455 ParamName);
456
457 // NOTE: absl::beta_distribution is not guaranteed to be stable.
TEST(BetaDistributionTest,StabilityTest)458 TEST(BetaDistributionTest, StabilityTest) {
459 // absl::beta_distribution stability relies on the stability of
460 // absl::random_interna::RandU64ToDouble, std::exp, std::log, std::pow,
461 // and std::sqrt.
462 //
463 // This test also depends on the stability of std::frexp.
464 using testing::ElementsAre;
465 absl::random_internal::sequence_urbg urbg({
466 0xffff00000000e6c8ull, 0xffff0000000006c8ull, 0x800003766295CFA9ull,
467 0x11C819684E734A41ull, 0x832603766295CFA9ull, 0x7fbe76c8b4395800ull,
468 0xB3472DCA7B14A94Aull, 0x0003eb76f6f7f755ull, 0xFFCEA50FDB2F953Bull,
469 0x13CCA830EB61BD96ull, 0x0334FE1EAA0363CFull, 0x00035C904C70A239ull,
470 0x00009E0BCBAADE14ull, 0x0000000000622CA7ull, 0x4864f22c059bf29eull,
471 0x247856d8b862665cull, 0xe46e86e9a1337e10ull, 0xd8c8541f3519b133ull,
472 0xffe75b52c567b9e4ull, 0xfffff732e5709c5bull, 0xff1f7f0b983532acull,
473 0x1ec2e8986d2362caull, 0xC332DDEFBE6C5AA5ull, 0x6558218568AB9702ull,
474 0x2AEF7DAD5B6E2F84ull, 0x1521B62829076170ull, 0xECDD4775619F1510ull,
475 0x814c8e35fe9a961aull, 0x0c3cd59c9b638a02ull, 0xcb3bb6478a07715cull,
476 0x1224e62c978bbc7full, 0x671ef2cb04e81f6eull, 0x3c1cbd811eaf1808ull,
477 0x1bbc23cfa8fac721ull, 0xa4c2cda65e596a51ull, 0xb77216fad37adf91ull,
478 0x836d794457c08849ull, 0xe083df03475f49d7ull, 0xbc9feb512e6b0d6cull,
479 0xb12d74fdd718c8c5ull, 0x12ff09653bfbe4caull, 0x8dd03a105bc4ee7eull,
480 0x5738341045ba0d85ull, 0xf3fd722dc65ad09eull, 0xfa14fd21ea2a5705ull,
481 0xffe6ea4d6edb0c73ull, 0xD07E9EFE2BF11FB4ull, 0x95DBDA4DAE909198ull,
482 0xEAAD8E716B93D5A0ull, 0xD08ED1D0AFC725E0ull, 0x8E3C5B2F8E7594B7ull,
483 0x8FF6E2FBF2122B64ull, 0x8888B812900DF01Cull, 0x4FAD5EA0688FC31Cull,
484 0xD1CFF191B3A8C1ADull, 0x2F2F2218BE0E1777ull, 0xEA752DFE8B021FA1ull,
485 });
486
487 // Convert the real-valued result into a unit64 where we compare
488 // 5 (float) or 10 (double) decimal digits plus the base-2 exponent.
489 auto float_to_u64 = [](float d) {
490 int exp = 0;
491 auto f = std::frexp(d, &exp);
492 return (static_cast<uint64_t>(1e5 * f) * 10000) + std::abs(exp);
493 };
494 auto double_to_u64 = [](double d) {
495 int exp = 0;
496 auto f = std::frexp(d, &exp);
497 return (static_cast<uint64_t>(1e10 * f) * 10000) + std::abs(exp);
498 };
499
500 std::vector<uint64_t> output(20);
501 {
502 // Algorithm Joehnk (float)
503 absl::beta_distribution<float> dist(0.1f, 0.2f);
504 std::generate(std::begin(output), std::end(output),
505 [&] { return float_to_u64(dist(urbg)); });
506 EXPECT_EQ(44, urbg.invocations());
507 EXPECT_THAT(output, //
508 testing::ElementsAre(
509 998340000, 619030004, 500000001, 999990000, 996280000,
510 500000001, 844740004, 847210001, 999970000, 872320000,
511 585480007, 933280000, 869080042, 647670031, 528240004,
512 969980004, 626050008, 915930002, 833440033, 878040015));
513 }
514
515 urbg.reset();
516 {
517 // Algorithm Joehnk (double)
518 absl::beta_distribution<double> dist(0.1, 0.2);
519 std::generate(std::begin(output), std::end(output),
520 [&] { return double_to_u64(dist(urbg)); });
521 EXPECT_EQ(44, urbg.invocations());
522 EXPECT_THAT(
523 output, //
524 testing::ElementsAre(
525 99834713000000, 61903356870004, 50000000000001, 99999721170000,
526 99628374770000, 99999999990000, 84474397860004, 84721276240001,
527 99997407490000, 87232528120000, 58548364780007, 93328932910000,
528 86908237770042, 64767917930031, 52824581970004, 96998544140004,
529 62605946270008, 91593604380002, 83345031740033, 87804397230015));
530 }
531
532 urbg.reset();
533 {
534 // Algorithm Cheng 1
535 absl::beta_distribution<double> dist(0.9, 2.0);
536 std::generate(std::begin(output), std::end(output),
537 [&] { return double_to_u64(dist(urbg)); });
538 EXPECT_EQ(62, urbg.invocations());
539 EXPECT_THAT(
540 output, //
541 testing::ElementsAre(
542 62069004780001, 64433204450001, 53607416560000, 89644295430008,
543 61434586310019, 55172615890002, 62187161490000, 56433684810003,
544 80454622050005, 86418558710003, 92920514700001, 64645184680001,
545 58549183380000, 84881283650005, 71078728590002, 69949694970000,
546 73157461710001, 68592191300001, 70747623900000, 78584696930005));
547 }
548
549 urbg.reset();
550 {
551 // Algorithm Cheng 2
552 absl::beta_distribution<double> dist(1.5, 2.5);
553 std::generate(std::begin(output), std::end(output),
554 [&] { return double_to_u64(dist(urbg)); });
555 EXPECT_EQ(54, urbg.invocations());
556 EXPECT_THAT(
557 output, //
558 testing::ElementsAre(
559 75000029250001, 76751482860001, 53264575220000, 69193133650005,
560 78028324470013, 91573587560002, 59167523770000, 60658618560002,
561 80075870540000, 94141320460004, 63196592770003, 78883906300002,
562 96797992590001, 76907587800001, 56645167560000, 65408302280003,
563 53401156320001, 64731238570000, 83065573750001, 79788333820001));
564 }
565 }
566
567 // This is an implementation-specific test. If any part of the implementation
568 // changes, then it is likely that this test will change as well. Also, if
569 // dependencies of the distribution change, such as RandU64ToDouble, then this
570 // is also likely to change.
TEST(BetaDistributionTest,AlgorithmBounds)571 TEST(BetaDistributionTest, AlgorithmBounds) {
572 #if (defined(__i386__) || defined(_M_IX86)) && FLT_EVAL_METHOD != 0
573 // We're using an x87-compatible FPU, and intermediate operations are
574 // performed with 80-bit floats. This produces slightly different results from
575 // what we expect below.
576 GTEST_SKIP()
577 << "Skipping the test because we detected x87 floating-point semantics";
578 #endif
579
580 {
581 absl::random_internal::sequence_urbg urbg(
582 {0x7fbe76c8b4395800ull, 0x8000000000000000ull});
583 // u=0.499, v=0.5
584 absl::beta_distribution<double> dist(1e-4, 1e-4);
585 double a = dist(urbg);
586 EXPECT_EQ(a, 2.0202860861567108529e-09);
587 EXPECT_EQ(2, urbg.invocations());
588 }
589
590 // Test that both the float & double algorithms appropriately reject the
591 // initial draw.
592 {
593 // 1/alpha = 1/beta = 2.
594 absl::beta_distribution<float> dist(0.5, 0.5);
595
596 // first two outputs are close to 1.0 - epsilon,
597 // thus: (u ^ 2 + v ^ 2) > 1.0
598 absl::random_internal::sequence_urbg urbg(
599 {0xffff00000006e6c8ull, 0xffff00000007c7c8ull, 0x800003766295CFA9ull,
600 0x11C819684E734A41ull});
601 {
602 double y = absl::beta_distribution<double>(0.5, 0.5)(urbg);
603 EXPECT_EQ(4, urbg.invocations());
604 EXPECT_EQ(y, 0.9810668952633862) << y;
605 }
606
607 // ...and: log(u) * a ~= log(v) * b ~= -0.02
608 // thus z ~= -0.02 + log(1 + e(~0))
609 // ~= -0.02 + 0.69
610 // thus z > 0
611 urbg.reset();
612 {
613 float x = absl::beta_distribution<float>(0.5, 0.5)(urbg);
614 EXPECT_EQ(4, urbg.invocations());
615 EXPECT_NEAR(0.98106688261032104, x, 0.0000005) << x << "f";
616 }
617 }
618 }
619
620 } // namespace
621