1 // Copyright (c) the JPEG XL Project Authors. All rights reserved.
2 //
3 // Use of this source code is governed by a BSD-style
4 // license that can be found in the LICENSE file.
5
6 #include "lib/jxl/base/descriptive_statistics.h"
7
8 #include <stdio.h>
9 #include <stdlib.h>
10
11 #include <random>
12 #include <vector>
13
14 #include "gtest/gtest.h"
15 #include "lib/jxl/base/compiler_specific.h"
16 #include "lib/jxl/noise_distributions.h"
17
18 namespace jxl {
19 namespace {
20
21 // Assigns x to one of two streams so we can later test Assimilate.
22 template <typename Random>
NotifyEither(float x,Random * rng,Stats * JXL_RESTRICT stats1,Stats * JXL_RESTRICT stats2)23 void NotifyEither(float x, Random* rng, Stats* JXL_RESTRICT stats1,
24 Stats* JXL_RESTRICT stats2) {
25 if ((*rng)() & 128) {
26 stats1->Notify(x);
27 } else {
28 stats2->Notify(x);
29 }
30 }
31
TEST(StatsTest,TestGaussian)32 TEST(StatsTest, TestGaussian) {
33 Stats stats;
34 Stats stats1, stats2;
35 const float mean = 5.0f;
36 const float stddev = 4.0f;
37 NoiseGaussian noise(stddev);
38 std::mt19937 rng(129);
39 for (size_t i = 0; i < 1000 * 1000; ++i) {
40 const float x = noise(mean, &rng);
41 stats.Notify(x);
42 NotifyEither(x, &rng, &stats1, &stats2);
43 }
44 EXPECT_NEAR(mean, stats.Mean(), 0.01);
45 EXPECT_NEAR(stddev, stats.StandardDeviation(), 0.02);
46 EXPECT_NEAR(0.0, stats.Skewness(), 0.02);
47 EXPECT_NEAR(0.0, stats.Kurtosis() - 3, 0.02);
48 printf("%s\n", stats.ToString().c_str());
49
50 // Same results after merging both accumulators.
51 stats1.Assimilate(stats2);
52 EXPECT_NEAR(mean, stats1.Mean(), 0.01);
53 EXPECT_NEAR(stddev, stats1.StandardDeviation(), 0.02);
54 EXPECT_NEAR(0.0, stats1.Skewness(), 0.02);
55 EXPECT_NEAR(0.0, stats1.Kurtosis() - 3, 0.02);
56 }
57
TEST(StatsTest,TestUniform)58 TEST(StatsTest, TestUniform) {
59 Stats stats;
60 Stats stats1, stats2;
61 NoiseUniform noise(0, 256);
62 std::mt19937 rng(129), rng_split(65537);
63 for (size_t i = 0; i < 1000 * 1000; ++i) {
64 const float x = noise(0.0f, &rng);
65 stats.Notify(x);
66 NotifyEither(x, &rng_split, &stats1, &stats2);
67 }
68 EXPECT_NEAR(128.0, stats.Mean(), 0.05);
69 EXPECT_NEAR(0.0, stats.Min(), 0.01);
70 EXPECT_NEAR(256.0, stats.Max(), 0.01);
71 EXPECT_NEAR(70, stats.StandardDeviation(), 10);
72 // No outliers.
73 EXPECT_NEAR(-1.2, stats.Kurtosis() - 3, 0.1);
74 printf("%s\n", stats.ToString().c_str());
75
76 // Same results after merging both accumulators.
77 stats1.Assimilate(stats2);
78 EXPECT_NEAR(128.0, stats1.Mean(), 0.05);
79 EXPECT_NEAR(0.0, stats1.Min(), 0.01);
80 EXPECT_NEAR(256.0, stats1.Max(), 0.01);
81 EXPECT_NEAR(70, stats1.StandardDeviation(), 10);
82 }
83
TEST(StatsTest,CompareCentralMomentsAgainstTwoPass)84 TEST(StatsTest, CompareCentralMomentsAgainstTwoPass) {
85 // Vary seed so the thresholds are not specific to one distribution.
86 for (int rep = 0; rep < 200; ++rep) {
87 // Uniform avoids outliers.
88 NoiseUniform noise(0, 256);
89 std::mt19937 rng(129 + 13 * rep), rng_split(65537);
90
91 // Small count so bias (population vs sample) is visible.
92 const size_t kSamples = 20;
93
94 // First pass: compute mean
95 std::vector<float> samples;
96 samples.reserve(kSamples);
97 double sum = 0.0;
98 for (size_t i = 0; i < kSamples; ++i) {
99 const float x = noise(0.0f, &rng);
100 samples.push_back(x);
101 sum += x;
102 }
103 const double mean = sum / kSamples;
104
105 // Second pass: compute stats and moments
106 Stats stats;
107 Stats stats1, stats2;
108 double sum2 = 0.0;
109 double sum3 = 0.0;
110 double sum4 = 0.0;
111 for (const double x : samples) {
112 const double d = x - mean;
113 sum2 += d * d;
114 sum3 += d * d * d;
115 sum4 += d * d * d * d;
116
117 stats.Notify(x);
118 NotifyEither(x, &rng_split, &stats1, &stats2);
119 }
120 const double mu1 = mean;
121 const double mu2 = sum2 / kSamples;
122 const double mu3 = sum3 / kSamples;
123 const double mu4 = sum4 / kSamples;
124
125 // Raw central moments (note: Mu1 is zero by definition)
126 EXPECT_NEAR(mu1, stats.Mu1(), 1E-13);
127 EXPECT_NEAR(mu2, stats.Mu2(), 1E-11);
128 EXPECT_NEAR(mu3, stats.Mu3(), 1E-9);
129 EXPECT_NEAR(mu4, stats.Mu4(), 1E-6);
130
131 // Same results after merging both accumulators.
132 stats1.Assimilate(stats2);
133 EXPECT_NEAR(mu1, stats1.Mu1(), 1E-13);
134 EXPECT_NEAR(mu2, stats1.Mu2(), 1E-11);
135 EXPECT_NEAR(mu3, stats1.Mu3(), 1E-9);
136 EXPECT_NEAR(mu4, stats1.Mu4(), 1E-6);
137
138 const double sample_variance = mu2;
139 // Scaling factor for sampling bias
140 const double r = (kSamples - 1.0) / kSamples;
141 const double skewness = mu3 * pow(r / mu2, 1.5);
142 const double kurtosis = mu4 * pow(r / mu2, 2.0);
143
144 EXPECT_NEAR(sample_variance, stats.SampleVariance(),
145 sample_variance * 1E-12);
146 EXPECT_NEAR(skewness, stats.Skewness(), std::abs(skewness * 1E-11));
147 EXPECT_NEAR(kurtosis, stats.Kurtosis(), kurtosis * 1E-12);
148 }
149 }
150
151 } // namespace
152 } // namespace jxl
153