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

llvm-svn: 103910

1 files changed, 68 insertions, 8 deletions

diff --git a/libcxx/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp b/libcxx/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp
index 4bb9ea10096..2632ab5f3ae 100644
--- a/libcxx/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp
+++ b/libcxx/test/numerics/rand/rand.dis/rand.dist.bern/rand.dist.bern.bin/eval_param.pass.cpp
@@ -39,60 +39,120 @@ int main()
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
+        {
+            D::result_type v = d(g, p);
+            assert(0 <= v && v <= p.t());
+            u.push_back(v);
+        }
         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.t() * p.p();
         double x_var = x_mean*(1-p.p());
+        double x_skew = (1-2*p.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*p.p()*(1-p.p())) / 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::binomial_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(16, .75);
         P p(30, .03125);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
+        {
+            D::result_type v = d(g, p);
+            assert(0 <= v && v <= p.t());
+            u.push_back(v);
+        }
         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.t() * p.p();
         double x_var = x_mean*(1-p.p());
+        double x_skew = (1-2*p.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*p.p()*(1-p.p())) / 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::binomial_distribution<> D;
         typedef D::param_type P;
-        typedef std::minstd_rand G;
+        typedef std::mt19937 G;
         G g;
         D d(16, .75);
         P p(40, .25);
         const int N = 100000;
         std::vector<D::result_type> u;
         for (int i = 0; i < N; ++i)
-            u.push_back(d(g, p));
+        {
+            D::result_type v = d(g, p);
+            assert(0 <= v && v <= p.t());
+            u.push_back(v);
+        }
         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.t() * p.p();
         double x_var = x_mean*(1-p.p());
+        double x_skew = (1-2*p.p()) / std::sqrt(x_var);
+        double x_kurtosis = (1-6*p.p()*(1-p.p())) / 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.03);
+        assert(std::abs(kurtosis - x_kurtosis) / x_kurtosis < 0.01);
     }
 }