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/gaussian_distribution.h"
16
17 #include <algorithm>
18 #include <cmath>
19 #include <cstddef>
20 #include <ios>
21 #include <iterator>
22 #include <random>
23 #include <string>
24 #include <vector>
25
26 #include "gmock/gmock.h"
27 #include "gtest/gtest.h"
28 #include "absl/base/internal/raw_logging.h"
29 #include "absl/base/macros.h"
30 #include "absl/random/internal/chi_square.h"
31 #include "absl/random/internal/distribution_test_util.h"
32 #include "absl/random/internal/sequence_urbg.h"
33 #include "absl/random/random.h"
34 #include "absl/strings/str_cat.h"
35 #include "absl/strings/str_format.h"
36 #include "absl/strings/str_replace.h"
37 #include "absl/strings/strip.h"
38
39 namespace {
40
41 using absl::random_internal::kChiSquared;
42
43 template <typename RealType>
44 class GaussianDistributionInterfaceTest : public ::testing::Test {};
45
46 using RealTypes = ::testing::Types<float, double, long double>;
47 TYPED_TEST_CASE(GaussianDistributionInterfaceTest, RealTypes);
48
TYPED_TEST(GaussianDistributionInterfaceTest,SerializeTest)49 TYPED_TEST(GaussianDistributionInterfaceTest, SerializeTest) {
50 using param_type =
51 typename absl::gaussian_distribution<TypeParam>::param_type;
52
53 const TypeParam kParams[] = {
54 // Cases around 1.
55 1, //
56 std::nextafter(TypeParam(1), TypeParam(0)), // 1 - epsilon
57 std::nextafter(TypeParam(1), TypeParam(2)), // 1 + epsilon
58 // Arbitrary values.
59 TypeParam(1e-8), TypeParam(1e-4), TypeParam(2), TypeParam(1e4),
60 TypeParam(1e8), TypeParam(1e20), TypeParam(2.5),
61 // Boundary cases.
62 std::numeric_limits<TypeParam>::infinity(),
63 std::numeric_limits<TypeParam>::max(),
64 std::numeric_limits<TypeParam>::epsilon(),
65 std::nextafter(std::numeric_limits<TypeParam>::min(),
66 TypeParam(1)), // min + epsilon
67 std::numeric_limits<TypeParam>::min(), // smallest normal
68 // There are some errors dealing with denorms on apple platforms.
69 std::numeric_limits<TypeParam>::denorm_min(), // smallest denorm
70 std::numeric_limits<TypeParam>::min() / 2,
71 std::nextafter(std::numeric_limits<TypeParam>::min(),
72 TypeParam(0)), // denorm_max
73 };
74
75 constexpr int kCount = 1000;
76 absl::InsecureBitGen gen;
77
78 // Use a loop to generate the combinations of {+/-x, +/-y}, and assign x, y to
79 // all values in kParams,
80 for (const auto mod : {0, 1, 2, 3}) {
81 for (const auto x : kParams) {
82 if (!std::isfinite(x)) continue;
83 for (const auto y : kParams) {
84 const TypeParam mean = (mod & 0x1) ? -x : x;
85 const TypeParam stddev = (mod & 0x2) ? -y : y;
86 const param_type param(mean, stddev);
87
88 absl::gaussian_distribution<TypeParam> before(mean, stddev);
89 EXPECT_EQ(before.mean(), param.mean());
90 EXPECT_EQ(before.stddev(), param.stddev());
91
92 {
93 absl::gaussian_distribution<TypeParam> via_param(param);
94 EXPECT_EQ(via_param, before);
95 EXPECT_EQ(via_param.param(), before.param());
96 }
97
98 // Smoke test.
99 auto sample_min = before.max();
100 auto sample_max = before.min();
101 for (int i = 0; i < kCount; i++) {
102 auto sample = before(gen);
103 if (sample > sample_max) sample_max = sample;
104 if (sample < sample_min) sample_min = sample;
105 EXPECT_GE(sample, before.min()) << before;
106 EXPECT_LE(sample, before.max()) << before;
107 }
108 if (!std::is_same<TypeParam, long double>::value) {
109 ABSL_INTERNAL_LOG(
110 INFO, absl::StrFormat("Range{%f, %f}: %f, %f", mean, stddev,
111 sample_min, sample_max));
112 }
113
114 std::stringstream ss;
115 ss << before;
116
117 if (!std::isfinite(mean) || !std::isfinite(stddev)) {
118 // Streams do not parse inf/nan.
119 continue;
120 }
121
122 // Validate stream serialization.
123 absl::gaussian_distribution<TypeParam> after(-0.53f, 2.3456f);
124
125 EXPECT_NE(before.mean(), after.mean());
126 EXPECT_NE(before.stddev(), after.stddev());
127 EXPECT_NE(before.param(), after.param());
128 EXPECT_NE(before, after);
129
130 ss >> after;
131
132 #if defined(__powerpc64__) || defined(__PPC64__) || defined(__powerpc__) || \
133 defined(__ppc__) || defined(__PPC__) || defined(__EMSCRIPTEN__)
134 if (std::is_same<TypeParam, long double>::value) {
135 // Roundtripping floating point values requires sufficient precision
136 // to reconstruct the exact value. It turns out that long double
137 // has some errors doing this on ppc, particularly for values
138 // near {1.0 +/- epsilon}.
139 //
140 // Emscripten is even worse, implementing long double as a 128-bit
141 // type, but shipping with a strtold() that doesn't support that.
142 if (mean <= std::numeric_limits<double>::max() &&
143 mean >= std::numeric_limits<double>::lowest()) {
144 EXPECT_EQ(static_cast<double>(before.mean()),
145 static_cast<double>(after.mean()))
146 << ss.str();
147 }
148 if (stddev <= std::numeric_limits<double>::max() &&
149 stddev >= std::numeric_limits<double>::lowest()) {
150 EXPECT_EQ(static_cast<double>(before.stddev()),
151 static_cast<double>(after.stddev()))
152 << ss.str();
153 }
154 continue;
155 }
156 #endif
157
158 EXPECT_EQ(before.mean(), after.mean());
159 EXPECT_EQ(before.stddev(), after.stddev()) //
160 << ss.str() << " " //
161 << (ss.good() ? "good " : "") //
162 << (ss.bad() ? "bad " : "") //
163 << (ss.eof() ? "eof " : "") //
164 << (ss.fail() ? "fail " : "");
165 }
166 }
167 }
168 }
169
170 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
171
172 class GaussianModel {
173 public:
GaussianModel(double mean,double stddev)174 GaussianModel(double mean, double stddev) : mean_(mean), stddev_(stddev) {}
175
mean() const176 double mean() const { return mean_; }
variance() const177 double variance() const { return stddev() * stddev(); }
stddev() const178 double stddev() const { return stddev_; }
skew() const179 double skew() const { return 0; }
kurtosis() const180 double kurtosis() const { return 3.0; }
181
182 // The inverse CDF, or PercentPoint function.
InverseCDF(double p)183 double InverseCDF(double p) {
184 ABSL_ASSERT(p >= 0.0);
185 ABSL_ASSERT(p < 1.0);
186 return mean() + stddev() * -absl::random_internal::InverseNormalSurvival(p);
187 }
188
189 private:
190 const double mean_;
191 const double stddev_;
192 };
193
194 struct Param {
195 double mean;
196 double stddev;
197 double p_fail; // Z-Test probability of failure.
198 int trials; // Z-Test trials.
199 };
200
201 // GaussianDistributionTests implements a z-test for the gaussian
202 // distribution.
203 class GaussianDistributionTests : public testing::TestWithParam<Param>,
204 public GaussianModel {
205 public:
GaussianDistributionTests()206 GaussianDistributionTests()
207 : GaussianModel(GetParam().mean, GetParam().stddev) {}
208
209 // SingleZTest provides a basic z-squared test of the mean vs. expected
210 // mean for data generated by the poisson distribution.
211 template <typename D>
212 bool SingleZTest(const double p, const size_t samples);
213
214 // SingleChiSquaredTest provides a basic chi-squared test of the normal
215 // distribution.
216 template <typename D>
217 double SingleChiSquaredTest();
218
219 // We use a fixed bit generator for distribution accuracy tests. This allows
220 // these tests to be deterministic, while still testing the qualify of the
221 // implementation.
222 absl::random_internal::pcg64_2018_engine rng_{0x2B7E151628AED2A6};
223 };
224
225 template <typename D>
SingleZTest(const double p,const size_t samples)226 bool GaussianDistributionTests::SingleZTest(const double p,
227 const size_t samples) {
228 D dis(mean(), stddev());
229
230 std::vector<double> data;
231 data.reserve(samples);
232 for (size_t i = 0; i < samples; i++) {
233 const double x = dis(rng_);
234 data.push_back(x);
235 }
236
237 const double max_err = absl::random_internal::MaxErrorTolerance(p);
238 const auto m = absl::random_internal::ComputeDistributionMoments(data);
239 const double z = absl::random_internal::ZScore(mean(), m);
240 const bool pass = absl::random_internal::Near("z", z, 0.0, max_err);
241
242 // NOTE: Informational statistical test:
243 //
244 // Compute the Jarque-Bera test statistic given the excess skewness
245 // and kurtosis. The statistic is drawn from a chi-square(2) distribution.
246 // https://en.wikipedia.org/wiki/Jarque%E2%80%93Bera_test
247 //
248 // The null-hypothesis (normal distribution) is rejected when
249 // (p = 0.05 => jb > 5.99)
250 // (p = 0.01 => jb > 9.21)
251 // NOTE: JB has a large type-I error rate, so it will reject the
252 // null-hypothesis even when it is true more often than the z-test.
253 //
254 const double jb =
255 static_cast<double>(m.n) / 6.0 *
256 (std::pow(m.skewness, 2.0) + std::pow(m.kurtosis - 3.0, 2.0) / 4.0);
257
258 if (!pass || jb > 9.21) {
259 ABSL_INTERNAL_LOG(
260 INFO, absl::StrFormat("p=%f max_err=%f\n"
261 " mean=%f vs. %f\n"
262 " stddev=%f vs. %f\n"
263 " skewness=%f vs. %f\n"
264 " kurtosis=%f vs. %f\n"
265 " z=%f vs. 0\n"
266 " jb=%f vs. 9.21",
267 p, max_err, m.mean, mean(), std::sqrt(m.variance),
268 stddev(), m.skewness, skew(), m.kurtosis,
269 kurtosis(), z, jb));
270 }
271 return pass;
272 }
273
274 template <typename D>
SingleChiSquaredTest()275 double GaussianDistributionTests::SingleChiSquaredTest() {
276 const size_t kSamples = 10000;
277 const int kBuckets = 50;
278
279 // The InverseCDF is the percent point function of the
280 // distribution, and can be used to assign buckets
281 // roughly uniformly.
282 std::vector<double> cutoffs;
283 const double kInc = 1.0 / static_cast<double>(kBuckets);
284 for (double p = kInc; p < 1.0; p += kInc) {
285 cutoffs.push_back(InverseCDF(p));
286 }
287 if (cutoffs.back() != std::numeric_limits<double>::infinity()) {
288 cutoffs.push_back(std::numeric_limits<double>::infinity());
289 }
290
291 D dis(mean(), stddev());
292
293 std::vector<int32_t> counts(cutoffs.size(), 0);
294 for (int j = 0; j < kSamples; j++) {
295 const double x = dis(rng_);
296 auto it = std::upper_bound(cutoffs.begin(), cutoffs.end(), x);
297 counts[std::distance(cutoffs.begin(), it)]++;
298 }
299
300 // Null-hypothesis is that the distribution is a gaussian distribution
301 // with the provided mean and stddev (not estimated from the data).
302 const int dof = static_cast<int>(counts.size()) - 1;
303
304 // Our threshold for logging is 1-in-50.
305 const double threshold = absl::random_internal::ChiSquareValue(dof, 0.98);
306
307 const double expected =
308 static_cast<double>(kSamples) / static_cast<double>(counts.size());
309
310 double chi_square = absl::random_internal::ChiSquareWithExpected(
311 std::begin(counts), std::end(counts), expected);
312 double p = absl::random_internal::ChiSquarePValue(chi_square, dof);
313
314 // Log if the chi_square value is above the threshold.
315 if (chi_square > threshold) {
316 for (int i = 0; i < cutoffs.size(); i++) {
317 ABSL_INTERNAL_LOG(
318 INFO, absl::StrFormat("%d : (%f) = %d", i, cutoffs[i], counts[i]));
319 }
320
321 ABSL_INTERNAL_LOG(
322 INFO, absl::StrCat("mean=", mean(), " stddev=", stddev(), "\n", //
323 " expected ", expected, "\n", //
324 kChiSquared, " ", chi_square, " (", p, ")\n", //
325 kChiSquared, " @ 0.98 = ", threshold));
326 }
327 return p;
328 }
329
TEST_P(GaussianDistributionTests,ZTest)330 TEST_P(GaussianDistributionTests, ZTest) {
331 // TODO(absl-team): Run these tests against std::normal_distribution<double>
332 // to validate outcomes are similar.
333 const size_t kSamples = 10000;
334 const auto& param = GetParam();
335 const int expected_failures =
336 std::max(1, static_cast<int>(std::ceil(param.trials * param.p_fail)));
337 const double p = absl::random_internal::RequiredSuccessProbability(
338 param.p_fail, param.trials);
339
340 int failures = 0;
341 for (int i = 0; i < param.trials; i++) {
342 failures +=
343 SingleZTest<absl::gaussian_distribution<double>>(p, kSamples) ? 0 : 1;
344 }
345 EXPECT_LE(failures, expected_failures);
346 }
347
TEST_P(GaussianDistributionTests,ChiSquaredTest)348 TEST_P(GaussianDistributionTests, ChiSquaredTest) {
349 const int kTrials = 20;
350 int failures = 0;
351
352 for (int i = 0; i < kTrials; i++) {
353 double p_value =
354 SingleChiSquaredTest<absl::gaussian_distribution<double>>();
355 if (p_value < 0.0025) { // 1/400
356 failures++;
357 }
358 }
359 // There is a 0.05% chance of producing at least one failure, so raise the
360 // failure threshold high enough to allow for a flake rate of less than one in
361 // 10,000.
362 EXPECT_LE(failures, 4);
363 }
364
GenParams()365 std::vector<Param> GenParams() {
366 return {
367 // Mean around 0.
368 Param{0.0, 1.0, 0.01, 100},
369 Param{0.0, 1e2, 0.01, 100},
370 Param{0.0, 1e4, 0.01, 100},
371 Param{0.0, 1e8, 0.01, 100},
372 Param{0.0, 1e16, 0.01, 100},
373 Param{0.0, 1e-3, 0.01, 100},
374 Param{0.0, 1e-5, 0.01, 100},
375 Param{0.0, 1e-9, 0.01, 100},
376 Param{0.0, 1e-17, 0.01, 100},
377
378 // Mean around 1.
379 Param{1.0, 1.0, 0.01, 100},
380 Param{1.0, 1e2, 0.01, 100},
381 Param{1.0, 1e-2, 0.01, 100},
382
383 // Mean around 100 / -100
384 Param{1e2, 1.0, 0.01, 100},
385 Param{-1e2, 1.0, 0.01, 100},
386 Param{1e2, 1e6, 0.01, 100},
387 Param{-1e2, 1e6, 0.01, 100},
388
389 // More extreme
390 Param{1e4, 1e4, 0.01, 100},
391 Param{1e8, 1e4, 0.01, 100},
392 Param{1e12, 1e4, 0.01, 100},
393 };
394 }
395
ParamName(const::testing::TestParamInfo<Param> & info)396 std::string ParamName(const ::testing::TestParamInfo<Param>& info) {
397 const auto& p = info.param;
398 std::string name = absl::StrCat("mean_", absl::SixDigits(p.mean), "__stddev_",
399 absl::SixDigits(p.stddev));
400 return absl::StrReplaceAll(name, {{"+", "_"}, {"-", "_"}, {".", "_"}});
401 }
402
403 INSTANTIATE_TEST_SUITE_P(All, GaussianDistributionTests,
404 ::testing::ValuesIn(GenParams()), ParamName);
405
406 // NOTE: absl::gaussian_distribution is not guaranteed to be stable.
TEST(GaussianDistributionTest,StabilityTest)407 TEST(GaussianDistributionTest, StabilityTest) {
408 // absl::gaussian_distribution stability relies on the underlying zignor
409 // data, absl::random_interna::RandU64ToDouble, std::exp, std::log, and
410 // std::abs.
411 absl::random_internal::sequence_urbg urbg(
412 {0x0003eb76f6f7f755ull, 0xFFCEA50FDB2F953Bull, 0xC332DDEFBE6C5AA5ull,
413 0x6558218568AB9702ull, 0x2AEF7DAD5B6E2F84ull, 0x1521B62829076170ull,
414 0xECDD4775619F1510ull, 0x13CCA830EB61BD96ull, 0x0334FE1EAA0363CFull,
415 0xB5735C904C70A239ull, 0xD59E9E0BCBAADE14ull, 0xEECC86BC60622CA7ull});
416
417 std::vector<int> output(11);
418
419 {
420 absl::gaussian_distribution<double> dist;
421 std::generate(std::begin(output), std::end(output),
422 [&] { return static_cast<int>(10000000.0 * dist(urbg)); });
423
424 EXPECT_EQ(13, urbg.invocations());
425 EXPECT_THAT(output, //
426 testing::ElementsAre(1494, 25518841, 9991550, 1351856,
427 -20373238, 3456682, 333530, -6804981,
428 -15279580, -16459654, 1494));
429 }
430
431 urbg.reset();
432 {
433 absl::gaussian_distribution<float> dist;
434 std::generate(std::begin(output), std::end(output),
435 [&] { return static_cast<int>(1000000.0f * dist(urbg)); });
436
437 EXPECT_EQ(13, urbg.invocations());
438 EXPECT_THAT(
439 output, //
440 testing::ElementsAre(149, 2551884, 999155, 135185, -2037323, 345668,
441 33353, -680498, -1527958, -1645965, 149));
442 }
443 }
444
445 // This is an implementation-specific test. If any part of the implementation
446 // changes, then it is likely that this test will change as well.
447 // Also, if dependencies of the distribution change, such as RandU64ToDouble,
448 // then this is also likely to change.
TEST(GaussianDistributionTest,AlgorithmBounds)449 TEST(GaussianDistributionTest, AlgorithmBounds) {
450 absl::gaussian_distribution<double> dist;
451
452 // In ~95% of cases, a single value is used to generate the output.
453 // for all inputs where |x| < 0.750461021389 this should be the case.
454 //
455 // The exact constraints are based on the ziggurat tables, and any
456 // changes to the ziggurat tables may require adjusting these bounds.
457 //
458 // for i in range(0, len(X)-1):
459 // print i, X[i+1]/X[i], (X[i+1]/X[i] > 0.984375)
460 //
461 // 0.125 <= |values| <= 0.75
462 const uint64_t kValues[] = {
463 0x1000000000000100ull, 0x2000000000000100ull, 0x3000000000000100ull,
464 0x4000000000000100ull, 0x5000000000000100ull, 0x6000000000000100ull,
465 // negative values
466 0x9000000000000100ull, 0xa000000000000100ull, 0xb000000000000100ull,
467 0xc000000000000100ull, 0xd000000000000100ull, 0xe000000000000100ull};
468
469 // 0.875 <= |values| <= 0.984375
470 const uint64_t kExtraValues[] = {
471 0x7000000000000100ull, 0x7800000000000100ull, //
472 0x7c00000000000100ull, 0x7e00000000000100ull, //
473 // negative values
474 0xf000000000000100ull, 0xf800000000000100ull, //
475 0xfc00000000000100ull, 0xfe00000000000100ull};
476
477 auto make_box = [](uint64_t v, uint64_t box) {
478 return (v & 0xffffffffffffff80ull) | box;
479 };
480
481 // The box is the lower 7 bits of the value. When the box == 0, then
482 // the algorithm uses an escape hatch to select the result for large
483 // outputs.
484 for (uint64_t box = 0; box < 0x7f; box++) {
485 for (const uint64_t v : kValues) {
486 // Extra values are added to the sequence to attempt to avoid
487 // infinite loops from rejection sampling on bugs/errors.
488 absl::random_internal::sequence_urbg urbg(
489 {make_box(v, box), 0x0003eb76f6f7f755ull, 0x5FCEA50FDB2F953Bull});
490
491 auto a = dist(urbg);
492 EXPECT_EQ(1, urbg.invocations()) << box << " " << std::hex << v;
493 if (v & 0x8000000000000000ull) {
494 EXPECT_LT(a, 0.0) << box << " " << std::hex << v;
495 } else {
496 EXPECT_GT(a, 0.0) << box << " " << std::hex << v;
497 }
498 }
499 if (box > 10 && box < 100) {
500 // The center boxes use the fast algorithm for more
501 // than 98.4375% of values.
502 for (const uint64_t v : kExtraValues) {
503 absl::random_internal::sequence_urbg urbg(
504 {make_box(v, box), 0x0003eb76f6f7f755ull, 0x5FCEA50FDB2F953Bull});
505
506 auto a = dist(urbg);
507 EXPECT_EQ(1, urbg.invocations()) << box << " " << std::hex << v;
508 if (v & 0x8000000000000000ull) {
509 EXPECT_LT(a, 0.0) << box << " " << std::hex << v;
510 } else {
511 EXPECT_GT(a, 0.0) << box << " " << std::hex << v;
512 }
513 }
514 }
515 }
516
517 // When the box == 0, the fallback algorithm uses a ratio of uniforms,
518 // which consumes 2 additional values from the urbg.
519 // Fallback also requires that the initial value be > 0.9271586026096681.
520 auto make_fallback = [](uint64_t v) { return (v & 0xffffffffffffff80ull); };
521
522 double tail[2];
523 {
524 // 0.9375
525 absl::random_internal::sequence_urbg urbg(
526 {make_fallback(0x7800000000000000ull), 0x13CCA830EB61BD96ull,
527 0x00000076f6f7f755ull});
528 tail[0] = dist(urbg);
529 EXPECT_EQ(3, urbg.invocations());
530 EXPECT_GT(tail[0], 0);
531 }
532 {
533 // -0.9375
534 absl::random_internal::sequence_urbg urbg(
535 {make_fallback(0xf800000000000000ull), 0x13CCA830EB61BD96ull,
536 0x00000076f6f7f755ull});
537 tail[1] = dist(urbg);
538 EXPECT_EQ(3, urbg.invocations());
539 EXPECT_LT(tail[1], 0);
540 }
541 EXPECT_EQ(tail[0], -tail[1]);
542 EXPECT_EQ(418610, static_cast<int64_t>(tail[0] * 100000.0));
543
544 // When the box != 0, the fallback algorithm computes a wedge function.
545 // Depending on the box, the threshold for varies as high as
546 // 0.991522480228.
547 {
548 // 0.9921875, 0.875
549 absl::random_internal::sequence_urbg urbg(
550 {make_box(0x7f00000000000000ull, 120), 0xe000000000000001ull,
551 0x13CCA830EB61BD96ull});
552 tail[0] = dist(urbg);
553 EXPECT_EQ(2, urbg.invocations());
554 EXPECT_GT(tail[0], 0);
555 }
556 {
557 // -0.9921875, 0.875
558 absl::random_internal::sequence_urbg urbg(
559 {make_box(0xff00000000000000ull, 120), 0xe000000000000001ull,
560 0x13CCA830EB61BD96ull});
561 tail[1] = dist(urbg);
562 EXPECT_EQ(2, urbg.invocations());
563 EXPECT_LT(tail[1], 0);
564 }
565 EXPECT_EQ(tail[0], -tail[1]);
566 EXPECT_EQ(61948, static_cast<int64_t>(tail[0] * 100000.0));
567
568 // Fallback rejected, try again.
569 {
570 // -0.9921875, 0.0625
571 absl::random_internal::sequence_urbg urbg(
572 {make_box(0xff00000000000000ull, 120), 0x1000000000000001,
573 make_box(0x1000000000000100ull, 50), 0x13CCA830EB61BD96ull});
574 dist(urbg);
575 EXPECT_EQ(3, urbg.invocations());
576 }
577 }
578
579 } // namespace
580