| // Copyright 2017 The Abseil Authors. |
| // |
| // Licensed under the Apache License, Version 2.0 (the "License"); |
| // you may not use this file except in compliance with the License. |
| // You may obtain a copy of the License at |
| // |
| // https://www.apache.org/licenses/LICENSE-2.0 |
| // |
| // Unless required by applicable law or agreed to in writing, software |
| // distributed under the License is distributed on an "AS IS" BASIS, |
| // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| // See the License for the specific language governing permissions and |
| // limitations under the License. |
| |
| #include "absl/random/beta_distribution.h" |
| |
| #include <algorithm> |
| #include <cfloat> |
| #include <cstddef> |
| #include <cstdint> |
| #include <iterator> |
| #include <random> |
| #include <sstream> |
| #include <string> |
| #include <type_traits> |
| #include <unordered_map> |
| #include <vector> |
| |
| #include "gmock/gmock.h" |
| #include "gtest/gtest.h" |
| #include "absl/base/internal/raw_logging.h" |
| #include "absl/numeric/internal/representation.h" |
| #include "absl/random/internal/chi_square.h" |
| #include "absl/random/internal/distribution_test_util.h" |
| #include "absl/random/internal/pcg_engine.h" |
| #include "absl/random/internal/sequence_urbg.h" |
| #include "absl/random/random.h" |
| #include "absl/strings/str_cat.h" |
| #include "absl/strings/str_format.h" |
| #include "absl/strings/str_replace.h" |
| #include "absl/strings/strip.h" |
| |
| namespace { |
| |
| template <typename IntType> |
| class BetaDistributionInterfaceTest : public ::testing::Test {}; |
| |
| // double-double arithmetic is not supported well by either GCC or Clang; see |
| // https://gcc.gnu.org/bugzilla/show_bug.cgi?id=99048, |
| // https://bugs.llvm.org/show_bug.cgi?id=49131, and |
| // https://bugs.llvm.org/show_bug.cgi?id=49132. Don't bother running these tests |
| // with double doubles until compiler support is better. |
| using RealTypes = |
| std::conditional<absl::numeric_internal::IsDoubleDouble(), |
| ::testing::Types<float, double>, |
| ::testing::Types<float, double, long double>>::type; |
| TYPED_TEST_SUITE(BetaDistributionInterfaceTest, RealTypes); |
| |
| TYPED_TEST(BetaDistributionInterfaceTest, SerializeTest) { |
| // The threshold for whether std::exp(1/a) is finite. |
| const TypeParam kSmallA = |
| 1.0f / std::log((std::numeric_limits<TypeParam>::max)()); |
| // The threshold for whether a * std::log(a) is finite. |
| const TypeParam kLargeA = |
| std::exp(std::log((std::numeric_limits<TypeParam>::max)()) - |
| std::log(std::log((std::numeric_limits<TypeParam>::max)()))); |
| using param_type = typename absl::beta_distribution<TypeParam>::param_type; |
| |
| constexpr int kCount = 1000; |
| absl::InsecureBitGen gen; |
| const TypeParam kValues[] = { |
| TypeParam(1e-20), TypeParam(1e-12), TypeParam(1e-8), TypeParam(1e-4), |
| TypeParam(1e-3), TypeParam(0.1), TypeParam(0.25), |
| std::nextafter(TypeParam(0.5), TypeParam(0)), // 0.5 - epsilon |
| std::nextafter(TypeParam(0.5), TypeParam(1)), // 0.5 + epsilon |
| TypeParam(0.5), TypeParam(1.0), // |
| std::nextafter(TypeParam(1), TypeParam(0)), // 1 - epsilon |
| std::nextafter(TypeParam(1), TypeParam(2)), // 1 + epsilon |
| TypeParam(12.5), TypeParam(1e2), TypeParam(1e8), TypeParam(1e12), |
| TypeParam(1e20), // |
| kSmallA, // |
| std::nextafter(kSmallA, TypeParam(0)), // |
| std::nextafter(kSmallA, TypeParam(1)), // |
| kLargeA, // |
| std::nextafter(kLargeA, TypeParam(0)), // |
| std::nextafter(kLargeA, std::numeric_limits<TypeParam>::max()), |
| // Boundary cases. |
| std::numeric_limits<TypeParam>::max(), |
| std::numeric_limits<TypeParam>::epsilon(), |
| std::nextafter(std::numeric_limits<TypeParam>::min(), |
| TypeParam(1)), // min + epsilon |
| std::numeric_limits<TypeParam>::min(), // smallest normal |
| std::numeric_limits<TypeParam>::denorm_min(), // smallest denorm |
| std::numeric_limits<TypeParam>::min() / 2, // denorm |
| std::nextafter(std::numeric_limits<TypeParam>::min(), |
| TypeParam(0)), // denorm_max |
| }; |
| for (TypeParam alpha : kValues) { |
| for (TypeParam beta : kValues) { |
| ABSL_INTERNAL_LOG( |
| INFO, absl::StrFormat("Smoke test for Beta(%a, %a)", alpha, beta)); |
| |
| param_type param(alpha, beta); |
| absl::beta_distribution<TypeParam> before(alpha, beta); |
| EXPECT_EQ(before.alpha(), param.alpha()); |
| EXPECT_EQ(before.beta(), param.beta()); |
| |
| { |
| absl::beta_distribution<TypeParam> via_param(param); |
| EXPECT_EQ(via_param, before); |
| EXPECT_EQ(via_param.param(), before.param()); |
| } |
| |
| // Smoke test. |
| for (int i = 0; i < kCount; ++i) { |
| auto sample = before(gen); |
| EXPECT_TRUE(std::isfinite(sample)); |
| EXPECT_GE(sample, before.min()); |
| EXPECT_LE(sample, before.max()); |
| } |
| |
| // Validate stream serialization. |
| std::stringstream ss; |
| ss << before; |
| absl::beta_distribution<TypeParam> after(3.8f, 1.43f); |
| EXPECT_NE(before.alpha(), after.alpha()); |
| EXPECT_NE(before.beta(), after.beta()); |
| EXPECT_NE(before.param(), after.param()); |
| EXPECT_NE(before, after); |
| |
| ss >> after; |
| |
| EXPECT_EQ(before.alpha(), after.alpha()); |
| EXPECT_EQ(before.beta(), after.beta()); |
| EXPECT_EQ(before, after) // |
| << ss.str() << " " // |
| << (ss.good() ? "good " : "") // |
| << (ss.bad() ? "bad " : "") // |
| << (ss.eof() ? "eof " : "") // |
| << (ss.fail() ? "fail " : ""); |
| } |
| } |
| } |
| |
| TYPED_TEST(BetaDistributionInterfaceTest, DegenerateCases) { |
| // We use a fixed bit generator for distribution accuracy tests. This allows |
| // these tests to be deterministic, while still testing the qualify of the |
| // implementation. |
| absl::random_internal::pcg64_2018_engine rng(0x2B7E151628AED2A6); |
| |
| // Extreme cases when the params are abnormal. |
| constexpr int kCount = 1000; |
| const TypeParam kSmallValues[] = { |
| std::numeric_limits<TypeParam>::min(), |
| std::numeric_limits<TypeParam>::denorm_min(), |
| std::nextafter(std::numeric_limits<TypeParam>::min(), |
| TypeParam(0)), // denorm_max |
| std::numeric_limits<TypeParam>::epsilon(), |
| }; |
| const TypeParam kLargeValues[] = { |
| std::numeric_limits<TypeParam>::max() * static_cast<TypeParam>(0.9999), |
| std::numeric_limits<TypeParam>::max() - 1, |
| std::numeric_limits<TypeParam>::max(), |
| }; |
| { |
| // Small alpha and beta. |
| // Useful WolframAlpha plots: |
| // * plot InverseBetaRegularized[x, 0.0001, 0.0001] from 0.495 to 0.505 |
| // * Beta[1.0, 0.0000001, 0.0000001] |
| // * Beta[0.9999, 0.0000001, 0.0000001] |
| for (TypeParam alpha : kSmallValues) { |
| for (TypeParam beta : kSmallValues) { |
| int zeros = 0; |
| int ones = 0; |
| absl::beta_distribution<TypeParam> d(alpha, beta); |
| for (int i = 0; i < kCount; ++i) { |
| TypeParam x = d(rng); |
| if (x == 0.0) { |
| zeros++; |
| } else if (x == 1.0) { |
| ones++; |
| } |
| } |
| EXPECT_EQ(ones + zeros, kCount); |
| if (alpha == beta) { |
| EXPECT_NE(ones, 0); |
| EXPECT_NE(zeros, 0); |
| } |
| } |
| } |
| } |
| { |
| // Small alpha, large beta. |
| // Useful WolframAlpha plots: |
| // * plot InverseBetaRegularized[x, 0.0001, 10000] from 0.995 to 1 |
| // * Beta[0, 0.0000001, 1000000] |
| // * Beta[0.001, 0.0000001, 1000000] |
| // * Beta[1, 0.0000001, 1000000] |
| for (TypeParam alpha : kSmallValues) { |
| for (TypeParam beta : kLargeValues) { |
| absl::beta_distribution<TypeParam> d(alpha, beta); |
| for (int i = 0; i < kCount; ++i) { |
| EXPECT_EQ(d(rng), 0.0); |
| } |
| } |
| } |
| } |
| { |
| // Large alpha, small beta. |
| // Useful WolframAlpha plots: |
| // * plot InverseBetaRegularized[x, 10000, 0.0001] from 0 to 0.001 |
| // * Beta[0.99, 1000000, 0.0000001] |
| // * Beta[1, 1000000, 0.0000001] |
| for (TypeParam alpha : kLargeValues) { |
| for (TypeParam beta : kSmallValues) { |
| absl::beta_distribution<TypeParam> d(alpha, beta); |
| for (int i = 0; i < kCount; ++i) { |
| EXPECT_EQ(d(rng), 1.0); |
| } |
| } |
| } |
| } |
| { |
| // Large alpha and beta. |
| absl::beta_distribution<TypeParam> d(std::numeric_limits<TypeParam>::max(), |
| std::numeric_limits<TypeParam>::max()); |
| for (int i = 0; i < kCount; ++i) { |
| EXPECT_EQ(d(rng), 0.5); |
| } |
| } |
| { |
| // Large alpha and beta but unequal. |
| absl::beta_distribution<TypeParam> d( |
| std::numeric_limits<TypeParam>::max(), |
| std::numeric_limits<TypeParam>::max() * 0.9999); |
| for (int i = 0; i < kCount; ++i) { |
| TypeParam x = d(rng); |
| EXPECT_NE(x, 0.5f); |
| EXPECT_FLOAT_EQ(x, 0.500025f); |
| } |
| } |
| } |
| |
| class BetaDistributionModel { |
| public: |
| explicit BetaDistributionModel(::testing::tuple<double, double> p) |
| : alpha_(::testing::get<0>(p)), beta_(::testing::get<1>(p)) {} |
| |
| double Mean() const { return alpha_ / (alpha_ + beta_); } |
| |
| double Variance() const { |
| return alpha_ * beta_ / (alpha_ + beta_ + 1) / (alpha_ + beta_) / |
| (alpha_ + beta_); |
| } |
| |
| double Kurtosis() const { |
| return 3 + 6 * |
| ((alpha_ - beta_) * (alpha_ - beta_) * (alpha_ + beta_ + 1) - |
| alpha_ * beta_ * (2 + alpha_ + beta_)) / |
| alpha_ / beta_ / (alpha_ + beta_ + 2) / (alpha_ + beta_ + 3); |
| } |
| |
| protected: |
| const double alpha_; |
| const double beta_; |
| }; |
| |
| class BetaDistributionTest |
| : public ::testing::TestWithParam<::testing::tuple<double, double>>, |
| public BetaDistributionModel { |
| public: |
| BetaDistributionTest() : BetaDistributionModel(GetParam()) {} |
| |
| protected: |
| template <class D> |
| bool SingleZTestOnMeanAndVariance(double p, size_t samples); |
| |
| template <class D> |
| bool SingleChiSquaredTest(double p, size_t samples, size_t buckets); |
| |
| absl::InsecureBitGen rng_; |
| }; |
| |
| template <class D> |
| bool BetaDistributionTest::SingleZTestOnMeanAndVariance(double p, |
| size_t samples) { |
| D dis(alpha_, beta_); |
| |
| std::vector<double> data; |
| data.reserve(samples); |
| for (size_t i = 0; i < samples; i++) { |
| const double variate = dis(rng_); |
| EXPECT_FALSE(std::isnan(variate)); |
| // Note that equality is allowed on both sides. |
| EXPECT_GE(variate, 0.0); |
| EXPECT_LE(variate, 1.0); |
| data.push_back(variate); |
| } |
| |
| // We validate that the sample mean and sample variance are indeed from a |
| // Beta distribution with the given shape parameters. |
| const auto m = absl::random_internal::ComputeDistributionMoments(data); |
| |
| // The variance of the sample mean is variance / n. |
| const double mean_stddev = std::sqrt(Variance() / static_cast<double>(m.n)); |
| |
| // The variance of the sample variance is (approximately): |
| // (kurtosis - 1) * variance^2 / n |
| const double variance_stddev = std::sqrt( |
| (Kurtosis() - 1) * Variance() * Variance() / static_cast<double>(m.n)); |
| // z score for the sample variance. |
| const double z_variance = (m.variance - Variance()) / variance_stddev; |
| |
| const double max_err = absl::random_internal::MaxErrorTolerance(p); |
| const double z_mean = absl::random_internal::ZScore(Mean(), m); |
| const bool pass = |
| absl::random_internal::Near("z", z_mean, 0.0, max_err) && |
| absl::random_internal::Near("z_variance", z_variance, 0.0, max_err); |
| if (!pass) { |
| ABSL_INTERNAL_LOG( |
| INFO, |
| absl::StrFormat( |
| "Beta(%f, %f), " |
| "mean: sample %f, expect %f, which is %f stddevs away, " |
| "variance: sample %f, expect %f, which is %f stddevs away.", |
| alpha_, beta_, m.mean, Mean(), |
| std::abs(m.mean - Mean()) / mean_stddev, m.variance, Variance(), |
| std::abs(m.variance - Variance()) / variance_stddev)); |
| } |
| return pass; |
| } |
| |
| template <class D> |
| bool BetaDistributionTest::SingleChiSquaredTest(double p, size_t samples, |
| size_t buckets) { |
| constexpr double kErr = 1e-7; |
| std::vector<double> cutoffs, expected; |
| const double bucket_width = 1.0 / static_cast<double>(buckets); |
| int i = 1; |
| int unmerged_buckets = 0; |
| for (; i < buckets; ++i) { |
| const double p = bucket_width * static_cast<double>(i); |
| const double boundary = |
| absl::random_internal::BetaIncompleteInv(alpha_, beta_, p); |
| // The intention is to add `boundary` to the list of `cutoffs`. It becomes |
| // problematic, however, when the boundary values are not monotone, due to |
| // numerical issues when computing the inverse regularized incomplete |
| // Beta function. In these cases, we merge that bucket with its previous |
| // neighbor and merge their expected counts. |
| if ((cutoffs.empty() && boundary < kErr) || |
| (!cutoffs.empty() && boundary <= cutoffs.back())) { |
| unmerged_buckets++; |
| continue; |
| } |
| if (boundary >= 1.0 - 1e-10) { |
| break; |
| } |
| cutoffs.push_back(boundary); |
| expected.push_back(static_cast<double>(1 + unmerged_buckets) * |
| bucket_width * static_cast<double>(samples)); |
| unmerged_buckets = 0; |
| } |
| cutoffs.push_back(std::numeric_limits<double>::infinity()); |
| // Merge all remaining buckets. |
| expected.push_back(static_cast<double>(buckets - i + 1) * bucket_width * |
| static_cast<double>(samples)); |
| // Make sure that we don't merge all the buckets, making this test |
| // meaningless. |
| EXPECT_GE(cutoffs.size(), 3) << alpha_ << ", " << beta_; |
| |
| D dis(alpha_, beta_); |
| |
| std::vector<int32_t> counts(cutoffs.size(), 0); |
| for (int i = 0; i < samples; i++) { |
| const double x = dis(rng_); |
| auto it = std::upper_bound(cutoffs.begin(), cutoffs.end(), x); |
| counts[std::distance(cutoffs.begin(), it)]++; |
| } |
| |
| // Null-hypothesis is that the distribution is beta distributed with the |
| // provided alpha, beta params (not estimated from the data). |
| const int dof = cutoffs.size() - 1; |
| |
| const double chi_square = absl::random_internal::ChiSquare( |
| counts.begin(), counts.end(), expected.begin(), expected.end()); |
| const bool pass = |
| (absl::random_internal::ChiSquarePValue(chi_square, dof) >= p); |
| if (!pass) { |
| for (int i = 0; i < cutoffs.size(); i++) { |
| ABSL_INTERNAL_LOG( |
| INFO, absl::StrFormat("cutoff[%d] = %f, actual count %d, expected %d", |
| i, cutoffs[i], counts[i], |
| static_cast<int>(expected[i]))); |
| } |
| |
| ABSL_INTERNAL_LOG( |
| INFO, absl::StrFormat( |
| "Beta(%f, %f) %s %f, p = %f", alpha_, beta_, |
| absl::random_internal::kChiSquared, chi_square, |
| absl::random_internal::ChiSquarePValue(chi_square, dof))); |
| } |
| return pass; |
| } |
| |
| TEST_P(BetaDistributionTest, TestSampleStatistics) { |
| static constexpr int kRuns = 20; |
| static constexpr double kPFail = 0.02; |
| const double p = |
| absl::random_internal::RequiredSuccessProbability(kPFail, kRuns); |
| static constexpr int kSampleCount = 10000; |
| static constexpr int kBucketCount = 100; |
| int failed = 0; |
| for (int i = 0; i < kRuns; ++i) { |
| if (!SingleZTestOnMeanAndVariance<absl::beta_distribution<double>>( |
| p, kSampleCount)) { |
| failed++; |
| } |
| if (!SingleChiSquaredTest<absl::beta_distribution<double>>( |
| 0.005, kSampleCount, kBucketCount)) { |
| failed++; |
| } |
| } |
| // Set so that the test is not flaky at --runs_per_test=10000 |
| EXPECT_LE(failed, 5); |
| } |
| |
| std::string ParamName( |
| const ::testing::TestParamInfo<::testing::tuple<double, double>>& info) { |
| std::string name = absl::StrCat("alpha_", ::testing::get<0>(info.param), |
| "__beta_", ::testing::get<1>(info.param)); |
| return absl::StrReplaceAll(name, {{"+", "_"}, {"-", "_"}, {".", "_"}}); |
| } |
| |
| INSTANTIATE_TEST_SUITE_P( |
| TestSampleStatisticsCombinations, BetaDistributionTest, |
| ::testing::Combine(::testing::Values(0.1, 0.2, 0.9, 1.1, 2.5, 10.0, 123.4), |
| ::testing::Values(0.1, 0.2, 0.9, 1.1, 2.5, 10.0, 123.4)), |
| ParamName); |
| |
| INSTANTIATE_TEST_SUITE_P( |
| TestSampleStatistics_SelectedPairs, BetaDistributionTest, |
| ::testing::Values(std::make_pair(0.5, 1000), std::make_pair(1000, 0.5), |
| std::make_pair(900, 1000), std::make_pair(10000, 20000), |
| std::make_pair(4e5, 2e7), std::make_pair(1e7, 1e5)), |
| ParamName); |
| |
| // NOTE: absl::beta_distribution is not guaranteed to be stable. |
| TEST(BetaDistributionTest, StabilityTest) { |
| // absl::beta_distribution stability relies on the stability of |
| // absl::random_interna::RandU64ToDouble, std::exp, std::log, std::pow, |
| // and std::sqrt. |
| // |
| // This test also depends on the stability of std::frexp. |
| using testing::ElementsAre; |
| absl::random_internal::sequence_urbg urbg({ |
| 0xffff00000000e6c8ull, 0xffff0000000006c8ull, 0x800003766295CFA9ull, |
| 0x11C819684E734A41ull, 0x832603766295CFA9ull, 0x7fbe76c8b4395800ull, |
| 0xB3472DCA7B14A94Aull, 0x0003eb76f6f7f755ull, 0xFFCEA50FDB2F953Bull, |
| 0x13CCA830EB61BD96ull, 0x0334FE1EAA0363CFull, 0x00035C904C70A239ull, |
| 0x00009E0BCBAADE14ull, 0x0000000000622CA7ull, 0x4864f22c059bf29eull, |
| 0x247856d8b862665cull, 0xe46e86e9a1337e10ull, 0xd8c8541f3519b133ull, |
| 0xffe75b52c567b9e4ull, 0xfffff732e5709c5bull, 0xff1f7f0b983532acull, |
| 0x1ec2e8986d2362caull, 0xC332DDEFBE6C5AA5ull, 0x6558218568AB9702ull, |
| 0x2AEF7DAD5B6E2F84ull, 0x1521B62829076170ull, 0xECDD4775619F1510ull, |
| 0x814c8e35fe9a961aull, 0x0c3cd59c9b638a02ull, 0xcb3bb6478a07715cull, |
| 0x1224e62c978bbc7full, 0x671ef2cb04e81f6eull, 0x3c1cbd811eaf1808ull, |
| 0x1bbc23cfa8fac721ull, 0xa4c2cda65e596a51ull, 0xb77216fad37adf91ull, |
| 0x836d794457c08849ull, 0xe083df03475f49d7ull, 0xbc9feb512e6b0d6cull, |
| 0xb12d74fdd718c8c5ull, 0x12ff09653bfbe4caull, 0x8dd03a105bc4ee7eull, |
| 0x5738341045ba0d85ull, 0xf3fd722dc65ad09eull, 0xfa14fd21ea2a5705ull, |
| 0xffe6ea4d6edb0c73ull, 0xD07E9EFE2BF11FB4ull, 0x95DBDA4DAE909198ull, |
| 0xEAAD8E716B93D5A0ull, 0xD08ED1D0AFC725E0ull, 0x8E3C5B2F8E7594B7ull, |
| 0x8FF6E2FBF2122B64ull, 0x8888B812900DF01Cull, 0x4FAD5EA0688FC31Cull, |
| 0xD1CFF191B3A8C1ADull, 0x2F2F2218BE0E1777ull, 0xEA752DFE8B021FA1ull, |
| }); |
| |
| // Convert the real-valued result into a unit64 where we compare |
| // 5 (float) or 10 (double) decimal digits plus the base-2 exponent. |
| auto float_to_u64 = [](float d) { |
| int exp = 0; |
| auto f = std::frexp(d, &exp); |
| return (static_cast<uint64_t>(1e5 * f) * 10000) + std::abs(exp); |
| }; |
| auto double_to_u64 = [](double d) { |
| int exp = 0; |
| auto f = std::frexp(d, &exp); |
| return (static_cast<uint64_t>(1e10 * f) * 10000) + std::abs(exp); |
| }; |
| |
| std::vector<uint64_t> output(20); |
| { |
| // Algorithm Joehnk (float) |
| absl::beta_distribution<float> dist(0.1f, 0.2f); |
| std::generate(std::begin(output), std::end(output), |
| [&] { return float_to_u64(dist(urbg)); }); |
| EXPECT_EQ(44, urbg.invocations()); |
| EXPECT_THAT(output, // |
| testing::ElementsAre( |
| 998340000, 619030004, 500000001, 999990000, 996280000, |
| 500000001, 844740004, 847210001, 999970000, 872320000, |
| 585480007, 933280000, 869080042, 647670031, 528240004, |
| 969980004, 626050008, 915930002, 833440033, 878040015)); |
| } |
| |
| urbg.reset(); |
| { |
| // Algorithm Joehnk (double) |
| absl::beta_distribution<double> dist(0.1, 0.2); |
| std::generate(std::begin(output), std::end(output), |
| [&] { return double_to_u64(dist(urbg)); }); |
| EXPECT_EQ(44, urbg.invocations()); |
| EXPECT_THAT( |
| output, // |
| testing::ElementsAre( |
| 99834713000000, 61903356870004, 50000000000001, 99999721170000, |
| 99628374770000, 99999999990000, 84474397860004, 84721276240001, |
| 99997407490000, 87232528120000, 58548364780007, 93328932910000, |
| 86908237770042, 64767917930031, 52824581970004, 96998544140004, |
| 62605946270008, 91593604380002, 83345031740033, 87804397230015)); |
| } |
| |
| urbg.reset(); |
| { |
| // Algorithm Cheng 1 |
| absl::beta_distribution<double> dist(0.9, 2.0); |
| std::generate(std::begin(output), std::end(output), |
| [&] { return double_to_u64(dist(urbg)); }); |
| EXPECT_EQ(62, urbg.invocations()); |
| EXPECT_THAT( |
| output, // |
| testing::ElementsAre( |
| 62069004780001, 64433204450001, 53607416560000, 89644295430008, |
| 61434586310019, 55172615890002, 62187161490000, 56433684810003, |
| 80454622050005, 86418558710003, 92920514700001, 64645184680001, |
| 58549183380000, 84881283650005, 71078728590002, 69949694970000, |
| 73157461710001, 68592191300001, 70747623900000, 78584696930005)); |
| } |
| |
| urbg.reset(); |
| { |
| // Algorithm Cheng 2 |
| absl::beta_distribution<double> dist(1.5, 2.5); |
| std::generate(std::begin(output), std::end(output), |
| [&] { return double_to_u64(dist(urbg)); }); |
| EXPECT_EQ(54, urbg.invocations()); |
| EXPECT_THAT( |
| output, // |
| testing::ElementsAre( |
| 75000029250001, 76751482860001, 53264575220000, 69193133650005, |
| 78028324470013, 91573587560002, 59167523770000, 60658618560002, |
| 80075870540000, 94141320460004, 63196592770003, 78883906300002, |
| 96797992590001, 76907587800001, 56645167560000, 65408302280003, |
| 53401156320001, 64731238570000, 83065573750001, 79788333820001)); |
| } |
| } |
| |
| // This is an implementation-specific test. If any part of the implementation |
| // changes, then it is likely that this test will change as well. Also, if |
| // dependencies of the distribution change, such as RandU64ToDouble, then this |
| // is also likely to change. |
| TEST(BetaDistributionTest, AlgorithmBounds) { |
| #if (defined(__i386__) || defined(_M_IX86)) && FLT_EVAL_METHOD != 0 |
| // We're using an x87-compatible FPU, and intermediate operations are |
| // performed with 80-bit floats. This produces slightly different results from |
| // what we expect below. |
| GTEST_SKIP() |
| << "Skipping the test because we detected x87 floating-point semantics"; |
| #endif |
| |
| { |
| absl::random_internal::sequence_urbg urbg( |
| {0x7fbe76c8b4395800ull, 0x8000000000000000ull}); |
| // u=0.499, v=0.5 |
| absl::beta_distribution<double> dist(1e-4, 1e-4); |
| double a = dist(urbg); |
| EXPECT_EQ(a, 2.0202860861567108529e-09); |
| EXPECT_EQ(2, urbg.invocations()); |
| } |
| |
| // Test that both the float & double algorithms appropriately reject the |
| // initial draw. |
| { |
| // 1/alpha = 1/beta = 2. |
| absl::beta_distribution<float> dist(0.5, 0.5); |
| |
| // first two outputs are close to 1.0 - epsilon, |
| // thus: (u ^ 2 + v ^ 2) > 1.0 |
| absl::random_internal::sequence_urbg urbg( |
| {0xffff00000006e6c8ull, 0xffff00000007c7c8ull, 0x800003766295CFA9ull, |
| 0x11C819684E734A41ull}); |
| { |
| double y = absl::beta_distribution<double>(0.5, 0.5)(urbg); |
| EXPECT_EQ(4, urbg.invocations()); |
| EXPECT_EQ(y, 0.9810668952633862) << y; |
| } |
| |
| // ...and: log(u) * a ~= log(v) * b ~= -0.02 |
| // thus z ~= -0.02 + log(1 + e(~0)) |
| // ~= -0.02 + 0.69 |
| // thus z > 0 |
| urbg.reset(); |
| { |
| float x = absl::beta_distribution<float>(0.5, 0.5)(urbg); |
| EXPECT_EQ(4, urbg.invocations()); |
| EXPECT_NEAR(0.98106688261032104, x, 0.0000005) << x << "f"; |
| } |
| } |
| } |
| |
| } // namespace |
