From ca233a7ea0d614fcb2ba0273ffe08ce1dae3727b Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sun, 13 Oct 2019 21:17:14 -0400 Subject: [PATCH] Runs test [CI SKIP] --- doc/graphs/expected_runs_above_threshold.svg | 2 + doc/graphs/runs_test_statistic.svg | 2 + doc/statistics/runs_test.qbk | 140 +++++++++++++++++-- include/boost/math/statistics/runs_test.hpp | 45 ++++-- test/test_runs_test.cpp | 25 ++++ 5 files changed, 192 insertions(+), 22 deletions(-) create mode 100644 doc/graphs/expected_runs_above_threshold.svg create mode 100644 doc/graphs/runs_test_statistic.svg diff --git a/doc/graphs/expected_runs_above_threshold.svg b/doc/graphs/expected_runs_above_threshold.svg new file mode 100644 index 000000000..a9c64db6b --- /dev/null +++ b/doc/graphs/expected_runs_above_threshold.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/graphs/runs_test_statistic.svg b/doc/graphs/runs_test_statistic.svg new file mode 100644 index 000000000..ad492f119 --- /dev/null +++ b/doc/graphs/runs_test_statistic.svg @@ -0,0 +1,2 @@ + + \ No newline at end of file diff --git a/doc/statistics/runs_test.qbk b/doc/statistics/runs_test.qbk index 42d927898..3d76c8ca9 100644 --- a/doc/statistics/runs_test.qbk +++ b/doc/statistics/runs_test.qbk @@ -15,17 +15,35 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) namespace boost::math::statistics { template -std::pair runs_above_threshold(RandomAccessContainer const & v, typename RandomAccessContainer::value_type threshold); +std::pair runs_above_and_below_threshold(RandomAccessContainer const & v, typename RandomAccessContainer::value_type threshold); template -std::pair runs_above_median(RandomAccessContainer const & v); +std::pair runs_above_and_below_median(RandomAccessContainer const & v); }}} ``` [heading Background] -The "runs above median test" tries to determine if a sequence is random by looking at the number of runs which exceed the median value of the sequence. +The "runs above and below median test" tries to determine if a sequence is random by looking at the number of consecutive values which exceed (or are below) the median of the sequence. +For example, if we are given data {5, 2, 0, 4, 7, 9, 10, 6, 1, 8, 3}, first we find the median (5), and transform the vector into a list of + and -'s, depending on whether the value is greater or less than the median. +(Values equal to the median we simply ignore-this is a convention we have inherited from the wonderful `randtests` package in R.) +Hence our data vector is transformed into + +{-,-,-,+,+,+,-,+,-} + +which is 5 runs, with /n/[sub a] = 5 values above and /n/[sub b] = 5 values below the median. +[@https://www.itl.nist.gov/div898/handbook/eda/section3/eda35d.htm NIST] tells us the expected number of runs and their variance: + +[$../graphs/expected_runs_above_threshold.svg] + +from which we derive the test statistic + +[$../graphs/runs_test_statistic.svg] + +whose distribution we approximate as normal to extract a /p/-value. + + An example usage is as follows: @@ -35,27 +53,123 @@ An example usage is as follows: #include std::random_device rd; -std::mt19937 gen{rd()}; -std::normal_distribution dis{0,1}; -std::vector v(1024); -for (auto & x : v) { - x = dis(gen); -} - -auto [t, p] = boost::math::statistics::runs_above_median(v); +std::vector v{5, 2, 0, 4, 7, 9, 10, 6, 1, 8, 3}; +auto [t, p] = boost::math::statistics::runs_above_and_below_median(v); +// t = -0.670820393249936919; p = 0.502334954360502017; ``` - -The test statistic is the first element of the pair, and the /p/-value is the second element. +We see that we have about a 50% chance of seeing this number of runs if the null hypothesis of randomness is true, and hence the assumption of randomness seems reasonable. +As always, the test statistic is the first element of the pair, and the /p/-value is the second element. [heading Performance] There are two cases: Where the threshold (typically the median) has already been computed, and the case where the mean and sample variance must be computed on the fly. +Computing the median is fairly expensive (requiring a call to `boost::math::statistics::median`), and since the order of the original data must be preserved, it must allocate. +If you believe your data to come from a distribution where the means and median coincide, or if you've already computed the median in the course of some other analysis, then you can get away with a call to `runs_above_and_below_threshold` via + +``` +Real threshold = 5; +auto [t, p] = boost::math::statistics::runs_above_and_below_threshold(v, threshold); +``` + +The performance differences between these two cases are obvious: ``` ---------------------------------------------- Benchmark Time ---------------------------------------------- +BMRunsAboveAndBelowMedian/8 260 ns bytes_per_second=118.421M/s +BMRunsAboveAndBelowMedian/16 318 ns bytes_per_second=192.797M/s +BMRunsAboveAndBelowMedian/32 417 ns bytes_per_second=303.509M/s +BMRunsAboveAndBelowMedian/64 625 ns bytes_per_second=390.578M/s +BMRunsAboveAndBelowMedian/128 743 ns bytes_per_second=657.827M/s +BMRunsAboveAndBelowMedian/256 1308 ns bytes_per_second=767.181M/s +BMRunsAboveAndBelowMedian/512 1896 ns bytes_per_second=1034.31M/s +BMRunsAboveAndBelowMedian/1024 6582 ns bytes_per_second=594.458M/s +BMRunsAboveAndBelowMedian/2048 26067 ns bytes_per_second=300.001M/s +BMRunsAboveAndBelowMedian/4096 62023 ns bytes_per_second=252.125M/s +BMRunsAboveAndBelowMedian/8192 124976 ns bytes_per_second=250.256M/s +BMRunsAboveAndBelowMedian/16384 242171 ns bytes_per_second=258.29M/s +BMRunsAboveAndBelowMedian/32768 528683 ns bytes_per_second=236.714M/s +BMRunsAboveAndBelowMedian/65536 965354 ns bytes_per_second=259.185M/s +BMRunsAboveAndBelowMedian/131072 2514693 ns bytes_per_second=199.068M/s +BMRunsAboveAndBelowMedian/262144 4223084 ns bytes_per_second=237.058M/s +BMRunsAboveAndBelowMedian/524288 8638963 ns bytes_per_second=231.755M/s +BMRunsAboveAndBelowMedian/1048576 16215682 ns bytes_per_second=246.995M/s +BMRunsAboveAndBelowMedian/2097152 39180496 ns bytes_per_second=204.443M/s +BMRunsAboveAndBelowMedian/4194304 82495779 ns bytes_per_second=194.307M/s +BMRunsAboveAndBelowMedian/8388608 142675936 ns bytes_per_second=224.547M/s +BMRunsAboveAndBelowMedian/16777216 287638068 ns bytes_per_second=223.088M/s +BMRunsAboveAndBelowMedian_BigO 17.25 N +BMRunsAboveAndBelowMedian/8 191 ns bytes_per_second=320.129M/s +BMRunsAboveAndBelowMedian/16 233 ns bytes_per_second=523.526M/s +BMRunsAboveAndBelowMedian/32 334 ns bytes_per_second=730.8M/s +BMRunsAboveAndBelowMedian/64 456 ns bytes_per_second=1070.93M/s +BMRunsAboveAndBelowMedian/128 688 ns bytes_per_second=1.38789G/s +BMRunsAboveAndBelowMedian/256 1257 ns bytes_per_second=1.51807G/s +BMRunsAboveAndBelowMedian/512 2663 ns bytes_per_second=1.43406G/s +BMRunsAboveAndBelowMedian/1024 4100 ns bytes_per_second=1.86266G/s +BMRunsAboveAndBelowMedian/2048 23493 ns bytes_per_second=665.851M/s +BMRunsAboveAndBelowMedian/4096 57968 ns bytes_per_second=539.551M/s +BMRunsAboveAndBelowMedian/8192 142272 ns bytes_per_second=439.746M/s +BMRunsAboveAndBelowMedian/16384 260948 ns bytes_per_second=479.639M/s +BMRunsAboveAndBelowMedian/32768 551577 ns bytes_per_second=453.623M/s +BMRunsAboveAndBelowMedian/65536 1056583 ns bytes_per_second=473.654M/s +BMRunsAboveAndBelowMedian/131072 2123956 ns bytes_per_second=471.35M/s +BMRunsAboveAndBelowMedian/262144 5028745 ns bytes_per_second=398.111M/s +BMRunsAboveAndBelowMedian/524288 10399212 ns bytes_per_second=384.981M/s +BMRunsAboveAndBelowMedian/1048576 23089767 ns bytes_per_second=348.496M/s +BMRunsAboveAndBelowMedian/2097152 37626884 ns bytes_per_second=425.962M/s +BMRunsAboveAndBelowMedian/4194304 79281747 ns bytes_per_second=404.088M/s +BMRunsAboveAndBelowMedian/8388608 172055781 ns bytes_per_second=373.391M/s +BMRunsAboveAndBelowMedian/16777216 391377449 ns bytes_per_second=332.01M/s +BMRunsAboveAndBelowMedian_BigO 22.52 N +BMRunsAboveAndBelowThreshold/8 41.6 ns bytes_per_second=739.55M/s +BMRunsAboveAndBelowThreshold/16 58.4 ns bytes_per_second=1050.48M/s +BMRunsAboveAndBelowThreshold/32 66.5 ns bytes_per_second=1.79606G/s +BMRunsAboveAndBelowThreshold/64 115 ns bytes_per_second=2.0762G/s +BMRunsAboveAndBelowThreshold/128 198 ns bytes_per_second=2.41515G/s +BMRunsAboveAndBelowThreshold/256 365 ns bytes_per_second=2.61328G/s +BMRunsAboveAndBelowThreshold/512 720 ns bytes_per_second=2.65053G/s +BMRunsAboveAndBelowThreshold/1024 1424 ns bytes_per_second=2.68123G/s +BMRunsAboveAndBelowThreshold/2048 3009 ns bytes_per_second=2.5379G/s +BMRunsAboveAndBelowThreshold/4096 16748 ns bytes_per_second=933.699M/s +BMRunsAboveAndBelowThreshold/8192 40190 ns bytes_per_second=778.105M/s +BMRunsAboveAndBelowThreshold/16384 86500 ns bytes_per_second=723.067M/s +BMRunsAboveAndBelowThreshold/32768 176692 ns bytes_per_second=708.108M/s +BMRunsAboveAndBelowThreshold/65536 356863 ns bytes_per_second=701.198M/s +BMRunsAboveAndBelowThreshold/131072 714807 ns bytes_per_second=700.08M/s +BMRunsAboveAndBelowThreshold/262144 1429078 ns bytes_per_second=700.415M/s +BMRunsAboveAndBelowThreshold/524288 2877227 ns bytes_per_second=695.785M/s +BMRunsAboveAndBelowThreshold/1048576 5795662 ns bytes_per_second=691.222M/s +BMRunsAboveAndBelowThreshold/2097152 11562715 ns bytes_per_second=692.427M/s +BMRunsAboveAndBelowThreshold/4194304 23364846 ns bytes_per_second=686.464M/s +BMRunsAboveAndBelowThreshold/8388608 46442540 ns bytes_per_second=689.871M/s +BMRunsAboveAndBelowThreshold/16777216 92284501 ns bytes_per_second=694.006M/s +BMRunsAboveAndBelowThreshold_BigO 5.51 N +BMRunsAboveAndBelowThreshold/8 45.1 ns bytes_per_second=1.32169G/s +BMRunsAboveAndBelowThreshold/16 53.6 ns bytes_per_second=2.22712G/s +BMRunsAboveAndBelowThreshold/32 71.4 ns bytes_per_second=3.34079G/s +BMRunsAboveAndBelowThreshold/64 112 ns bytes_per_second=4.24946G/s +BMRunsAboveAndBelowThreshold/128 196 ns bytes_per_second=4.87317G/s +BMRunsAboveAndBelowThreshold/256 378 ns bytes_per_second=5.04476G/s +BMRunsAboveAndBelowThreshold/512 702 ns bytes_per_second=5.44134G/s +BMRunsAboveAndBelowThreshold/1024 1417 ns bytes_per_second=5.3865G/s +BMRunsAboveAndBelowThreshold/2048 3031 ns bytes_per_second=5.03872G/s +BMRunsAboveAndBelowThreshold/4096 16813 ns bytes_per_second=1.81669G/s +BMRunsAboveAndBelowThreshold/8192 41182 ns bytes_per_second=1.48565G/s +BMRunsAboveAndBelowThreshold/16384 86939 ns bytes_per_second=1.40536G/s +BMRunsAboveAndBelowThreshold/32768 177255 ns bytes_per_second=1.37892G/s +BMRunsAboveAndBelowThreshold/65536 356391 ns bytes_per_second=1.3713G/s +BMRunsAboveAndBelowThreshold/131072 718417 ns bytes_per_second=1.36057G/s +BMRunsAboveAndBelowThreshold/262144 1442288 ns bytes_per_second=1.35583G/s +BMRunsAboveAndBelowThreshold/524288 2942259 ns bytes_per_second=1.33217G/s +BMRunsAboveAndBelowThreshold/1048576 5870235 ns bytes_per_second=1.33244G/s +BMRunsAboveAndBelowThreshold/2097152 11743081 ns bytes_per_second=1.33192G/s +BMRunsAboveAndBelowThreshold/4194304 23521002 ns bytes_per_second=1.32976G/s +BMRunsAboveAndBelowThreshold/8388608 46917407 ns bytes_per_second=1.33339G/s +BMRunsAboveAndBelowThreshold/16777216 93823876 ns bytes_per_second=1.33305G/s +BMRunsAboveAndBelowThreshold_BigO 5.59 N 5.59 N ``` diff --git a/include/boost/math/statistics/runs_test.hpp b/include/boost/math/statistics/runs_test.hpp index 3e9846c12..dc1a6ecf6 100644 --- a/include/boost/math/statistics/runs_test.hpp +++ b/include/boost/math/statistics/runs_test.hpp @@ -10,6 +10,7 @@ #include #include +#include #include #include @@ -22,9 +23,9 @@ auto runs_above_and_below_threshold(RandomAccessContainer const & v, using Real = typename RandomAccessContainer::value_type; using std::sqrt; using std::abs; - if (v.size() <= 4) + if (v.size() <= 1) { - throw std::domain_error("At least 5 samples are required to get number of runs."); + throw std::domain_error("At least 2 samples are required to get number of runs."); } typedef boost::math::policies::policy< boost::math::policies::promote_float, @@ -34,16 +35,28 @@ auto runs_above_and_below_threshold(RandomAccessContainer const & v, decltype(v.size()) nabove = 0; decltype(v.size()) nbelow = 0; - bool run_up = (v[0] > threshold); + decltype(v.size()) imin = 0; + + // Take care of the case that v[0] == threshold: + while (imin < v.size() && v[imin] == threshold) { + ++imin; + } + + // Take care of the constant vector case: + if (imin == v.size()) { + return std::make_pair(std::numeric_limits::quiet_NaN(), Real(0)); + } + + bool run_up = (v[imin] > threshold); if (run_up) { ++nabove; } else { ++nbelow; } decltype(v.size()) runs = 1; - for (decltype(v.size()) i = 1; i < v.size(); ++i) { + for (decltype(v.size()) i = imin + 1; i < v.size(); ++i) { if (v[i] == threshold) { - // skip values precisely equal to threshold. + // skip values precisely equal to threshold (following R's randtests package) continue; } bool above = (v[i] > threshold); @@ -61,14 +74,28 @@ auto runs_above_and_below_threshold(RandomAccessContainer const & v, } } - // What should be done about a constant vector? Then n = 0: + // If you make n an int, the subtraction is gonna be bad in the variance: Real n = nabove + nbelow; Real expected_runs = Real(1) + Real(2*nabove*nbelow)/Real(n); - // What if this <= 0? I shudder to think of it. . . - Real var = 2*nabove*nbelow*(2*nabove*nbelow-n)/Real(n*n*(n-1)); + Real variance = 2*nabove*nbelow*(2*nabove*nbelow-n)/Real(n*n*(n-1)); - Real sd = sqrt(var); + // Bizarre, pathological limits: + if (variance == 0) + { + if (runs == expected_runs) + { + Real statistic = 0; + Real pvalue = 1; + return std::make_pair(statistic, pvalue); + } + else + { + return std::make_pair(std::numeric_limits::quiet_NaN(), Real(0)); + } + } + + Real sd = sqrt(variance); Real statistic = (runs - expected_runs)/sd; auto normal = boost::math::normal_distribution(0,1); diff --git a/test/test_runs_test.cpp b/test/test_runs_test.cpp index 450b7450e..94d5c5a2c 100644 --- a/test/test_runs_test.cpp +++ b/test/test_runs_test.cpp @@ -45,8 +45,33 @@ void test_agreement_with_r_randtests() CHECK_ULP_CLOSE(expected_pvalue, computed_pvalue, 3); } +void test_doc_example() +{ + std::vector v{5, 2, 0, 4, 7, 9, 10, 6, 1, 8, 3}; + double expected_statistic = -0.670820393249936919; + double expected_pvalue = 0.502334954360502017; + + auto [computed_statistic, computed_pvalue] = runs_above_and_below_median(v); + + CHECK_ULP_CLOSE(expected_statistic, computed_statistic, 3); + CHECK_ULP_CLOSE(expected_pvalue, computed_pvalue, 3); +} + +void test_constant_vector() +{ + std::vector v{5,5,5,5,5,5,5}; + auto [computed_statistic, computed_pvalue] = runs_above_and_below_median(v); + double expected_pvalue = 0; + CHECK_ULP_CLOSE(expected_pvalue, computed_pvalue, 3); + if (!std::isnan(computed_statistic)) { + std::cerr << "Computed statistic is not a nan!\n"; + } +} + int main() { + test_constant_vector(); test_agreement_with_r_randtests(); + test_doc_example(); return boost::math::test::report_errors(); }