boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Runs test [CI SKIP]

Browse Source
This commit is contained in:
Nick Thompson
2019-10-13 21:17:14 -04:00
parent c437f61911
commit ca233a7ea0
5 changed files with 192 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/graphs/expected_runs_above_threshold.svg Normal file
View File

File diff suppressed because one or more lines are too long
Side by Side

After

Width:  |  Height:  |  Size: 14 KiB

2
doc/graphs/runs_test_statistic.svg Normal file
View File

File diff suppressed because one or more lines are too long
Side by Side

After

Width:  |  Height:  |  Size: 5.1 KiB

140
doc/statistics/runs_test.qbk
View File

@@ -15,17 +15,35 @@ LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
namespace boost::math::statistics {
template<typename RandomAccessContainer>
std::pair<Real, Real> runs_above_threshold(RandomAccessContainer const & v, typename RandomAccessContainer::value_type threshold);
std::pair<Real, Real> runs_above_and_below_threshold(RandomAccessContainer const & v, typename RandomAccessContainer::value_type threshold);
template<typename RandomAccessContainer>
std::pair<Real, Real> runs_above_median(RandomAccessContainer const & v);
std::pair<Real, Real> 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 <boost/math/statistics/runs_test.hpp>
std::random_device rd;
std::mt19937 gen{rd()};
std::normal_distribution<double> dis{0,1};
std::vector<double> v(1024);
for (auto & x : v) {
x = dis(gen);
}
auto [t, p] = boost::math::statistics::runs_above_median(v);
std::vector<double> 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<float>/8 260 ns bytes_per_second=118.421M/s
BMRunsAboveAndBelowMedian<float>/16 318 ns bytes_per_second=192.797M/s
BMRunsAboveAndBelowMedian<float>/32 417 ns bytes_per_second=303.509M/s
BMRunsAboveAndBelowMedian<float>/64 625 ns bytes_per_second=390.578M/s
BMRunsAboveAndBelowMedian<float>/128 743 ns bytes_per_second=657.827M/s
BMRunsAboveAndBelowMedian<float>/256 1308 ns bytes_per_second=767.181M/s
BMRunsAboveAndBelowMedian<float>/512 1896 ns bytes_per_second=1034.31M/s
BMRunsAboveAndBelowMedian<float>/1024 6582 ns bytes_per_second=594.458M/s
BMRunsAboveAndBelowMedian<float>/2048 26067 ns bytes_per_second=300.001M/s
BMRunsAboveAndBelowMedian<float>/4096 62023 ns bytes_per_second=252.125M/s
BMRunsAboveAndBelowMedian<float>/8192 124976 ns bytes_per_second=250.256M/s
BMRunsAboveAndBelowMedian<float>/16384 242171 ns bytes_per_second=258.29M/s
BMRunsAboveAndBelowMedian<float>/32768 528683 ns bytes_per_second=236.714M/s
BMRunsAboveAndBelowMedian<float>/65536 965354 ns bytes_per_second=259.185M/s
BMRunsAboveAndBelowMedian<float>/131072 2514693 ns bytes_per_second=199.068M/s
BMRunsAboveAndBelowMedian<float>/262144 4223084 ns bytes_per_second=237.058M/s
BMRunsAboveAndBelowMedian<float>/524288 8638963 ns bytes_per_second=231.755M/s
BMRunsAboveAndBelowMedian<float>/1048576 16215682 ns bytes_per_second=246.995M/s
BMRunsAboveAndBelowMedian<float>/2097152 39180496 ns bytes_per_second=204.443M/s
BMRunsAboveAndBelowMedian<float>/4194304 82495779 ns bytes_per_second=194.307M/s
BMRunsAboveAndBelowMedian<float>/8388608 142675936 ns bytes_per_second=224.547M/s
BMRunsAboveAndBelowMedian<float>/16777216 287638068 ns bytes_per_second=223.088M/s
BMRunsAboveAndBelowMedian<float>_BigO 17.25 N
BMRunsAboveAndBelowMedian<double>/8 191 ns bytes_per_second=320.129M/s
BMRunsAboveAndBelowMedian<double>/16 233 ns bytes_per_second=523.526M/s
BMRunsAboveAndBelowMedian<double>/32 334 ns bytes_per_second=730.8M/s
BMRunsAboveAndBelowMedian<double>/64 456 ns bytes_per_second=1070.93M/s
BMRunsAboveAndBelowMedian<double>/128 688 ns bytes_per_second=1.38789G/s
BMRunsAboveAndBelowMedian<double>/256 1257 ns bytes_per_second=1.51807G/s
BMRunsAboveAndBelowMedian<double>/512 2663 ns bytes_per_second=1.43406G/s
BMRunsAboveAndBelowMedian<double>/1024 4100 ns bytes_per_second=1.86266G/s
BMRunsAboveAndBelowMedian<double>/2048 23493 ns bytes_per_second=665.851M/s
BMRunsAboveAndBelowMedian<double>/4096 57968 ns bytes_per_second=539.551M/s
BMRunsAboveAndBelowMedian<double>/8192 142272 ns bytes_per_second=439.746M/s
BMRunsAboveAndBelowMedian<double>/16384 260948 ns bytes_per_second=479.639M/s
BMRunsAboveAndBelowMedian<double>/32768 551577 ns bytes_per_second=453.623M/s
BMRunsAboveAndBelowMedian<double>/65536 1056583 ns bytes_per_second=473.654M/s
BMRunsAboveAndBelowMedian<double>/131072 2123956 ns bytes_per_second=471.35M/s
BMRunsAboveAndBelowMedian<double>/262144 5028745 ns bytes_per_second=398.111M/s
BMRunsAboveAndBelowMedian<double>/524288 10399212 ns bytes_per_second=384.981M/s
BMRunsAboveAndBelowMedian<double>/1048576 23089767 ns bytes_per_second=348.496M/s
BMRunsAboveAndBelowMedian<double>/2097152 37626884 ns bytes_per_second=425.962M/s
BMRunsAboveAndBelowMedian<double>/4194304 79281747 ns bytes_per_second=404.088M/s
BMRunsAboveAndBelowMedian<double>/8388608 172055781 ns bytes_per_second=373.391M/s
BMRunsAboveAndBelowMedian<double>/16777216 391377449 ns bytes_per_second=332.01M/s
BMRunsAboveAndBelowMedian<double>_BigO 22.52 N
BMRunsAboveAndBelowThreshold<float>/8 41.6 ns bytes_per_second=739.55M/s
BMRunsAboveAndBelowThreshold<float>/16 58.4 ns bytes_per_second=1050.48M/s
BMRunsAboveAndBelowThreshold<float>/32 66.5 ns bytes_per_second=1.79606G/s
BMRunsAboveAndBelowThreshold<float>/64 115 ns bytes_per_second=2.0762G/s
BMRunsAboveAndBelowThreshold<float>/128 198 ns bytes_per_second=2.41515G/s
BMRunsAboveAndBelowThreshold<float>/256 365 ns bytes_per_second=2.61328G/s
BMRunsAboveAndBelowThreshold<float>/512 720 ns bytes_per_second=2.65053G/s
BMRunsAboveAndBelowThreshold<float>/1024 1424 ns bytes_per_second=2.68123G/s
BMRunsAboveAndBelowThreshold<float>/2048 3009 ns bytes_per_second=2.5379G/s
BMRunsAboveAndBelowThreshold<float>/4096 16748 ns bytes_per_second=933.699M/s
BMRunsAboveAndBelowThreshold<float>/8192 40190 ns bytes_per_second=778.105M/s
BMRunsAboveAndBelowThreshold<float>/16384 86500 ns bytes_per_second=723.067M/s
BMRunsAboveAndBelowThreshold<float>/32768 176692 ns bytes_per_second=708.108M/s
BMRunsAboveAndBelowThreshold<float>/65536 356863 ns bytes_per_second=701.198M/s
BMRunsAboveAndBelowThreshold<float>/131072 714807 ns bytes_per_second=700.08M/s
BMRunsAboveAndBelowThreshold<float>/262144 1429078 ns bytes_per_second=700.415M/s
BMRunsAboveAndBelowThreshold<float>/524288 2877227 ns bytes_per_second=695.785M/s
BMRunsAboveAndBelowThreshold<float>/1048576 5795662 ns bytes_per_second=691.222M/s
BMRunsAboveAndBelowThreshold<float>/2097152 11562715 ns bytes_per_second=692.427M/s
BMRunsAboveAndBelowThreshold<float>/4194304 23364846 ns bytes_per_second=686.464M/s
BMRunsAboveAndBelowThreshold<float>/8388608 46442540 ns bytes_per_second=689.871M/s
BMRunsAboveAndBelowThreshold<float>/16777216 92284501 ns bytes_per_second=694.006M/s
BMRunsAboveAndBelowThreshold<float>_BigO 5.51 N
BMRunsAboveAndBelowThreshold<double>/8 45.1 ns bytes_per_second=1.32169G/s
BMRunsAboveAndBelowThreshold<double>/16 53.6 ns bytes_per_second=2.22712G/s
BMRunsAboveAndBelowThreshold<double>/32 71.4 ns bytes_per_second=3.34079G/s
BMRunsAboveAndBelowThreshold<double>/64 112 ns bytes_per_second=4.24946G/s
BMRunsAboveAndBelowThreshold<double>/128 196 ns bytes_per_second=4.87317G/s
BMRunsAboveAndBelowThreshold<double>/256 378 ns bytes_per_second=5.04476G/s
BMRunsAboveAndBelowThreshold<double>/512 702 ns bytes_per_second=5.44134G/s
BMRunsAboveAndBelowThreshold<double>/1024 1417 ns bytes_per_second=5.3865G/s
BMRunsAboveAndBelowThreshold<double>/2048 3031 ns bytes_per_second=5.03872G/s
BMRunsAboveAndBelowThreshold<double>/4096 16813 ns bytes_per_second=1.81669G/s
BMRunsAboveAndBelowThreshold<double>/8192 41182 ns bytes_per_second=1.48565G/s
BMRunsAboveAndBelowThreshold<double>/16384 86939 ns bytes_per_second=1.40536G/s
BMRunsAboveAndBelowThreshold<double>/32768 177255 ns bytes_per_second=1.37892G/s
BMRunsAboveAndBelowThreshold<double>/65536 356391 ns bytes_per_second=1.3713G/s
BMRunsAboveAndBelowThreshold<double>/131072 718417 ns bytes_per_second=1.36057G/s
BMRunsAboveAndBelowThreshold<double>/262144 1442288 ns bytes_per_second=1.35583G/s
BMRunsAboveAndBelowThreshold<double>/524288 2942259 ns bytes_per_second=1.33217G/s
BMRunsAboveAndBelowThreshold<double>/1048576 5870235 ns bytes_per_second=1.33244G/s
BMRunsAboveAndBelowThreshold<double>/2097152 11743081 ns bytes_per_second=1.33192G/s
BMRunsAboveAndBelowThreshold<double>/4194304 23521002 ns bytes_per_second=1.32976G/s
BMRunsAboveAndBelowThreshold<double>/8388608 46917407 ns bytes_per_second=1.33339G/s
BMRunsAboveAndBelowThreshold<double>/16777216 93823876 ns bytes_per_second=1.33305G/s
BMRunsAboveAndBelowThreshold<double>_BigO 5.59 N 5.59 N
```

45
include/boost/math/statistics/runs_test.hpp
View File

@@ -10,6 +10,7 @@
#include <cmath>
#include <algorithm>
#include <utility>
#include <boost/math/statistics/univariate_statistics.hpp>
#include <boost/math/distributions/normal.hpp>
@@ -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<false>,
@@ -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<Real>::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<Real>::quiet_NaN(), Real(0));
}
}
Real sd = sqrt(variance);
Real statistic = (runs - expected_runs)/sd;
auto normal = boost::math::normal_distribution<Real, no_promote_policy>(0,1);

25
test/test_runs_test.cpp
View File

@@ -45,8 +45,33 @@ void test_agreement_with_r_randtests()
CHECK_ULP_CLOSE(expected_pvalue, computed_pvalue, 3);
}
void test_doc_example()
{
std::vector<double> 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<double> 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();
}