beta_distribution_test.cc (23825B)
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/log/log.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 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 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 LOG(INFO) << absl::StreamFormat("Smoke test for Beta(%a, %a)", alpha, 111 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 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: 256 explicit BetaDistributionModel(::testing::tuple<double, double> p) 257 : alpha_(::testing::get<0>(p)), beta_(::testing::get<1>(p)) {} 258 259 double Mean() const { return alpha_ / (alpha_ + beta_); } 260 261 double Variance() const { 262 return alpha_ * beta_ / (alpha_ + beta_ + 1) / (alpha_ + beta_) / 263 (alpha_ + beta_); 264 } 265 266 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: 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> 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 LOG(INFO) << "Beta(" << alpha_ << ", " << beta_ << "), mean: sample " 331 << m.mean << ", expect " << Mean() << ", which is " 332 << std::abs(m.mean - Mean()) / mean_stddev 333 << " stddevs away, variance: sample " << m.variance << ", expect " 334 << Variance() << ", which is " 335 << std::abs(m.variance - Variance()) / variance_stddev 336 << " stddevs away."; 337 } 338 return pass; 339 } 340 341 template <class D> 342 bool BetaDistributionTest::SingleChiSquaredTest(double p, size_t samples, 343 size_t buckets) { 344 constexpr double kErr = 1e-7; 345 std::vector<double> cutoffs, expected; 346 const double bucket_width = 1.0 / static_cast<double>(buckets); 347 int i = 1; 348 int unmerged_buckets = 0; 349 for (; i < buckets; ++i) { 350 const double p = bucket_width * static_cast<double>(i); 351 const double boundary = 352 absl::random_internal::BetaIncompleteInv(alpha_, beta_, p); 353 // The intention is to add `boundary` to the list of `cutoffs`. It becomes 354 // problematic, however, when the boundary values are not monotone, due to 355 // numerical issues when computing the inverse regularized incomplete 356 // Beta function. In these cases, we merge that bucket with its previous 357 // neighbor and merge their expected counts. 358 if ((cutoffs.empty() && boundary < kErr) || 359 (!cutoffs.empty() && boundary <= cutoffs.back())) { 360 unmerged_buckets++; 361 continue; 362 } 363 if (boundary >= 1.0 - 1e-10) { 364 break; 365 } 366 cutoffs.push_back(boundary); 367 expected.push_back(static_cast<double>(1 + unmerged_buckets) * 368 bucket_width * static_cast<double>(samples)); 369 unmerged_buckets = 0; 370 } 371 cutoffs.push_back(std::numeric_limits<double>::infinity()); 372 // Merge all remaining buckets. 373 expected.push_back(static_cast<double>(buckets - i + 1) * bucket_width * 374 static_cast<double>(samples)); 375 // Make sure that we don't merge all the buckets, making this test 376 // meaningless. 377 EXPECT_GE(cutoffs.size(), 3) << alpha_ << ", " << beta_; 378 379 D dis(alpha_, beta_); 380 381 std::vector<int32_t> counts(cutoffs.size(), 0); 382 for (int i = 0; i < samples; i++) { 383 const double x = dis(rng_); 384 auto it = std::upper_bound(cutoffs.begin(), cutoffs.end(), x); 385 counts[std::distance(cutoffs.begin(), it)]++; 386 } 387 388 // Null-hypothesis is that the distribution is beta distributed with the 389 // provided alpha, beta params (not estimated from the data). 390 const int dof = cutoffs.size() - 1; 391 392 const double chi_square = absl::random_internal::ChiSquare( 393 counts.begin(), counts.end(), expected.begin(), expected.end()); 394 const bool pass = 395 (absl::random_internal::ChiSquarePValue(chi_square, dof) >= p); 396 if (!pass) { 397 for (size_t i = 0; i < cutoffs.size(); i++) { 398 LOG(INFO) << "cutoff[" << i << "] = " << cutoffs[i] << ", actual count " 399 << counts[i] << ", expected " << static_cast<int>(expected[i]); 400 } 401 402 LOG(INFO) << "Beta(" << alpha_ << ", " << beta_ << ") " 403 << absl::random_internal::kChiSquared << " " << chi_square 404 << ", p = " 405 << absl::random_internal::ChiSquarePValue(chi_square, dof); 406 } 407 return pass; 408 } 409 410 TEST_P(BetaDistributionTest, TestSampleStatistics) { 411 static constexpr int kRuns = 20; 412 static constexpr double kPFail = 0.02; 413 const double p = 414 absl::random_internal::RequiredSuccessProbability(kPFail, kRuns); 415 static constexpr int kSampleCount = 10000; 416 static constexpr int kBucketCount = 100; 417 int failed = 0; 418 for (int i = 0; i < kRuns; ++i) { 419 if (!SingleZTestOnMeanAndVariance<absl::beta_distribution<double>>( 420 p, kSampleCount)) { 421 failed++; 422 } 423 if (!SingleChiSquaredTest<absl::beta_distribution<double>>( 424 0.005, kSampleCount, kBucketCount)) { 425 failed++; 426 } 427 } 428 // Set so that the test is not flaky at --runs_per_test=10000 429 EXPECT_LE(failed, 5); 430 } 431 432 std::string ParamName( 433 const ::testing::TestParamInfo<::testing::tuple<double, double>>& info) { 434 std::string name = absl::StrCat("alpha_", ::testing::get<0>(info.param), 435 "__beta_", ::testing::get<1>(info.param)); 436 return absl::StrReplaceAll(name, {{"+", "_"}, {"-", "_"}, {".", "_"}}); 437 } 438 439 INSTANTIATE_TEST_SUITE_P( 440 TestSampleStatisticsCombinations, BetaDistributionTest, 441 ::testing::Combine(::testing::Values(0.1, 0.2, 0.9, 1.1, 2.5, 10.0, 123.4), 442 ::testing::Values(0.1, 0.2, 0.9, 1.1, 2.5, 10.0, 123.4)), 443 ParamName); 444 445 INSTANTIATE_TEST_SUITE_P( 446 TestSampleStatistics_SelectedPairs, BetaDistributionTest, 447 ::testing::Values(std::make_pair(0.5, 1000), std::make_pair(1000, 0.5), 448 std::make_pair(900, 1000), std::make_pair(10000, 20000), 449 std::make_pair(4e5, 2e7), std::make_pair(1e7, 1e5)), 450 ParamName); 451 452 // NOTE: absl::beta_distribution is not guaranteed to be stable. 453 TEST(BetaDistributionTest, StabilityTest) { 454 // absl::beta_distribution stability relies on the stability of 455 // absl::random_interna::RandU64ToDouble, std::exp, std::log, std::pow, 456 // and std::sqrt. 457 // 458 // This test also depends on the stability of std::frexp. 459 using testing::ElementsAre; 460 absl::random_internal::sequence_urbg urbg({ 461 0xffff00000000e6c8ull, 0xffff0000000006c8ull, 0x800003766295CFA9ull, 462 0x11C819684E734A41ull, 0x832603766295CFA9ull, 0x7fbe76c8b4395800ull, 463 0xB3472DCA7B14A94Aull, 0x0003eb76f6f7f755ull, 0xFFCEA50FDB2F953Bull, 464 0x13CCA830EB61BD96ull, 0x0334FE1EAA0363CFull, 0x00035C904C70A239ull, 465 0x00009E0BCBAADE14ull, 0x0000000000622CA7ull, 0x4864f22c059bf29eull, 466 0x247856d8b862665cull, 0xe46e86e9a1337e10ull, 0xd8c8541f3519b133ull, 467 0xffe75b52c567b9e4ull, 0xfffff732e5709c5bull, 0xff1f7f0b983532acull, 468 0x1ec2e8986d2362caull, 0xC332DDEFBE6C5AA5ull, 0x6558218568AB9702ull, 469 0x2AEF7DAD5B6E2F84ull, 0x1521B62829076170ull, 0xECDD4775619F1510ull, 470 0x814c8e35fe9a961aull, 0x0c3cd59c9b638a02ull, 0xcb3bb6478a07715cull, 471 0x1224e62c978bbc7full, 0x671ef2cb04e81f6eull, 0x3c1cbd811eaf1808ull, 472 0x1bbc23cfa8fac721ull, 0xa4c2cda65e596a51ull, 0xb77216fad37adf91ull, 473 0x836d794457c08849ull, 0xe083df03475f49d7ull, 0xbc9feb512e6b0d6cull, 474 0xb12d74fdd718c8c5ull, 0x12ff09653bfbe4caull, 0x8dd03a105bc4ee7eull, 475 0x5738341045ba0d85ull, 0xf3fd722dc65ad09eull, 0xfa14fd21ea2a5705ull, 476 0xffe6ea4d6edb0c73ull, 0xD07E9EFE2BF11FB4ull, 0x95DBDA4DAE909198ull, 477 0xEAAD8E716B93D5A0ull, 0xD08ED1D0AFC725E0ull, 0x8E3C5B2F8E7594B7ull, 478 0x8FF6E2FBF2122B64ull, 0x8888B812900DF01Cull, 0x4FAD5EA0688FC31Cull, 479 0xD1CFF191B3A8C1ADull, 0x2F2F2218BE0E1777ull, 0xEA752DFE8B021FA1ull, 480 }); 481 482 // Convert the real-valued result into a unit64 where we compare 483 // 5 (float) or 10 (double) decimal digits plus the base-2 exponent. 484 auto float_to_u64 = [](float d) { 485 int exp = 0; 486 auto f = std::frexp(d, &exp); 487 return (static_cast<uint64_t>(1e5 * f) * 10000) + std::abs(exp); 488 }; 489 auto double_to_u64 = [](double d) { 490 int exp = 0; 491 auto f = std::frexp(d, &exp); 492 return (static_cast<uint64_t>(1e10 * f) * 10000) + std::abs(exp); 493 }; 494 495 std::vector<uint64_t> output(20); 496 { 497 // Algorithm Joehnk (float) 498 absl::beta_distribution<float> dist(0.1f, 0.2f); 499 std::generate(std::begin(output), std::end(output), 500 [&] { return float_to_u64(dist(urbg)); }); 501 EXPECT_EQ(44, urbg.invocations()); 502 EXPECT_THAT(output, // 503 testing::ElementsAre( 504 998340000, 619030004, 500000001, 999990000, 996280000, 505 500000001, 844740004, 847210001, 999970000, 872320000, 506 585480007, 933280000, 869080042, 647670031, 528240004, 507 969980004, 626050008, 915930002, 833440033, 878040015)); 508 } 509 510 urbg.reset(); 511 { 512 // Algorithm Joehnk (double) 513 absl::beta_distribution<double> dist(0.1, 0.2); 514 std::generate(std::begin(output), std::end(output), 515 [&] { return double_to_u64(dist(urbg)); }); 516 EXPECT_EQ(44, urbg.invocations()); 517 EXPECT_THAT( 518 output, // 519 testing::ElementsAre( 520 99834713000000, 61903356870004, 50000000000001, 99999721170000, 521 99628374770000, 99999999990000, 84474397860004, 84721276240001, 522 99997407490000, 87232528120000, 58548364780007, 93328932910000, 523 86908237770042, 64767917930031, 52824581970004, 96998544140004, 524 62605946270008, 91593604380002, 83345031740033, 87804397230015)); 525 } 526 527 urbg.reset(); 528 { 529 // Algorithm Cheng 1 530 absl::beta_distribution<double> dist(0.9, 2.0); 531 std::generate(std::begin(output), std::end(output), 532 [&] { return double_to_u64(dist(urbg)); }); 533 EXPECT_EQ(62, urbg.invocations()); 534 EXPECT_THAT( 535 output, // 536 testing::ElementsAre( 537 62069004780001, 64433204450001, 53607416560000, 89644295430008, 538 61434586310019, 55172615890002, 62187161490000, 56433684810003, 539 80454622050005, 86418558710003, 92920514700001, 64645184680001, 540 58549183380000, 84881283650005, 71078728590002, 69949694970000, 541 73157461710001, 68592191300001, 70747623900000, 78584696930005)); 542 } 543 544 urbg.reset(); 545 { 546 // Algorithm Cheng 2 547 absl::beta_distribution<double> dist(1.5, 2.5); 548 std::generate(std::begin(output), std::end(output), 549 [&] { return double_to_u64(dist(urbg)); }); 550 EXPECT_EQ(54, urbg.invocations()); 551 EXPECT_THAT( 552 output, // 553 testing::ElementsAre( 554 75000029250001, 76751482860001, 53264575220000, 69193133650005, 555 78028324470013, 91573587560002, 59167523770000, 60658618560002, 556 80075870540000, 94141320460004, 63196592770003, 78883906300002, 557 96797992590001, 76907587800001, 56645167560000, 65408302280003, 558 53401156320001, 64731238570000, 83065573750001, 79788333820001)); 559 } 560 } 561 562 // This is an implementation-specific test. If any part of the implementation 563 // changes, then it is likely that this test will change as well. Also, if 564 // dependencies of the distribution change, such as RandU64ToDouble, then this 565 // is also likely to change. 566 TEST(BetaDistributionTest, AlgorithmBounds) { 567 #if (defined(__i386__) || defined(_M_IX86)) && FLT_EVAL_METHOD != 0 568 // We're using an x87-compatible FPU, and intermediate operations are 569 // performed with 80-bit floats. This produces slightly different results from 570 // what we expect below. 571 GTEST_SKIP() 572 << "Skipping the test because we detected x87 floating-point semantics"; 573 #endif 574 575 { 576 absl::random_internal::sequence_urbg urbg( 577 {0x7fbe76c8b4395800ull, 0x8000000000000000ull}); 578 // u=0.499, v=0.5 579 absl::beta_distribution<double> dist(1e-4, 1e-4); 580 double a = dist(urbg); 581 EXPECT_EQ(a, 2.0202860861567108529e-09); 582 EXPECT_EQ(2, urbg.invocations()); 583 } 584 585 // Test that both the float & double algorithms appropriately reject the 586 // initial draw. 587 { 588 // 1/alpha = 1/beta = 2. 589 absl::beta_distribution<float> dist(0.5, 0.5); 590 591 // first two outputs are close to 1.0 - epsilon, 592 // thus: (u ^ 2 + v ^ 2) > 1.0 593 absl::random_internal::sequence_urbg urbg( 594 {0xffff00000006e6c8ull, 0xffff00000007c7c8ull, 0x800003766295CFA9ull, 595 0x11C819684E734A41ull}); 596 { 597 double y = absl::beta_distribution<double>(0.5, 0.5)(urbg); 598 EXPECT_EQ(4, urbg.invocations()); 599 EXPECT_EQ(y, 0.9810668952633862) << y; 600 } 601 602 // ...and: log(u) * a ~= log(v) * b ~= -0.02 603 // thus z ~= -0.02 + log(1 + e(~0)) 604 // ~= -0.02 + 0.69 605 // thus z > 0 606 urbg.reset(); 607 { 608 float x = absl::beta_distribution<float>(0.5, 0.5)(urbg); 609 EXPECT_EQ(4, urbg.invocations()); 610 EXPECT_NEAR(0.98106688261032104, x, 0.0000005) << x << "f"; 611 } 612 } 613 } 614 615 } // namespace