Beefed up the tests for all of the distributions to include checks against the expected skewness and kurtosis

llvm-svn: 103910

Author: Howard Hinnant

libcxx/include/random

@@ -3434,9 +3434,7 @@ public:
     void param(const param_type& __p) {__p_ = __p;}
 
     result_type min() const {return 0;}
-    result_type max() const
-        {return -std::log(1-std::nextafter(result_type(1), result_type(-1))) /
-                 __p_.lambda();}
+    result_type max() const {return numeric_limits<result_type>::infinity();}
 
     friend bool operator==(const exponential_distribution& __x,
                            const exponential_distribution& __y)

libcxx/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval.pass.cpp

@@ -40,13 +40,29 @@ int main()
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.p();
         double x_var = d.p()*(1-d.p());
+        double x_skew = (1 - 2 * d.p())/std::sqrt(x_var);
+        double x_kurtosis = (6 * sqr(d.p()) - 6 * d.p() + 1)/x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::bernoulli_distribution D;
@@ -60,12 +76,28 @@ int main()
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = d.p();
         double x_var = d.p()*(1-d.p());
+        double x_skew = (1 - 2 * d.p())/std::sqrt(x_var);
+        double x_kurtosis = (6 * sqr(d.p()) - 6 * d.p() + 1)/x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }

libcxx/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bernoulli/eval_param.pass.cpp

@@ -42,13 +42,29 @@ int main()
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = p.p();
         double x_var = p.p()*(1-p.p());
+        double x_skew = (1 - 2 * p.p())/std::sqrt(x_var);
+        double x_kurtosis = (6 * sqr(p.p()) - 6 * p.p() + 1)/x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
     {
         typedef std::bernoulli_distribution D;
@@ -64,12 +80,28 @@ int main()
         double mean = std::accumulate(u.begin(), u.end(),
                                               double(0)) / u.size();
         double var = 0;
+        double skew = 0;
+        double kurtosis = 0;
         for (int i = 0; i < u.size(); ++i)
-            var += sqr(u[i] - mean);
+        {
+            double d = (u[i] - mean);
+            double d2 = sqr(d);
+            var += d2;
+            skew += d * d2;
+            kurtosis += d2 * d2;
+        }
         var /= u.size();
+        double dev = std::sqrt(var);
+        skew /= u.size() * dev * var;
+        kurtosis /= u.size() * var * var;
+        kurtosis -= 3;
         double x_mean = p.p();
         double x_var = p.p()*(1-p.p());
+        double x_skew = (1 - 2 * p.p())/std::sqrt(x_var);
+        double x_kurtosis = (6 * sqr(p.p()) - 6 * p.p() + 1)/x_var;
         assert(std::abs(mean - x_mean) / x_mean < 0.01);
         assert(std::abs(var - x_var) / x_var < 0.01);
+        assert(std::abs(skew - x_skew) / x_skew < 0.01);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }