math/statistics/runs_test.hpp

0001 /*
0002  * Copyright Nick Thompson, 2019
0003  * Use, modification and distribution are subject to the
0004  * Boost Software License, Version 1.0. (See accompanying file
0005  * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006  */
0007
0008 #ifndef BOOST_MATH_STATISTICS_RUNS_TEST_HPP
0009 #define BOOST_MATH_STATISTICS_RUNS_TEST_HPP
0010
0011 #include <cmath>
0012 #include <algorithm>
0013 #include <utility>
0014 #include <boost/math/statistics/univariate_statistics.hpp>
0015 #include <boost/math/distributions/normal.hpp>
0016
0017 namespace boost::math::statistics {
0018
0019 template<class RandomAccessContainer>
0020 auto runs_above_and_below_threshold(RandomAccessContainer const & v,
0021                           typename RandomAccessContainer::value_type threshold)
0022 {
0023     using Real = typename RandomAccessContainer::value_type;
0024     using std::sqrt;
0025     using std::abs;
0026     if (v.size() <= 1)
0027     {
0028         throw std::domain_error("At least 2 samples are required to get number of runs.");
0029     }
0030     typedef boost::math::policies::policy<
0031           boost::math::policies::promote_float<false>,
0032           boost::math::policies::promote_double<false> >
0033           no_promote_policy;
0034
0035     decltype(v.size()) nabove = 0;
0036     decltype(v.size()) nbelow = 0;
0037
0038     decltype(v.size()) imin = 0;
0039
0040     // Take care of the case that v[0] == threshold:
0041     while (imin < v.size() && v[imin] == threshold) {
0042         ++imin;
0043     }
0044
0045     // Take care of the constant vector case:
0046     if (imin == v.size()) {
0047         return std::make_pair(std::numeric_limits<Real>::quiet_NaN(), Real(0));
0048     }
0049
0050     bool run_up = (v[imin] > threshold);
0051     if (run_up) {
0052         ++nabove;
0053     } else {
0054         ++nbelow;
0055     }
0056     decltype(v.size()) runs = 1;
0057     for (decltype(v.size()) i = imin + 1; i < v.size(); ++i) {
0058       if (v[i] == threshold) {
0059         // skip values precisely equal to threshold (following R's randtests package)
0060         continue;
0061       }
0062       bool above = (v[i] > threshold);
0063       if (above) {
0064           ++nabove;
0065       } else {
0066           ++nbelow;
0067       }
0068       if (run_up == above) {
0069         continue;
0070       }
0071       else {
0072         run_up = above;
0073         runs++;
0074       }
0075     }
0076
0077     // If you make n an int, the subtraction is gonna be bad in the variance:
0078     Real n = nabove + nbelow;
0079
0080     Real expected_runs = Real(1) + Real(2*nabove*nbelow)/Real(n);
0081     Real variance = 2*nabove*nbelow*(2*nabove*nbelow-n)/Real(n*n*(n-1));
0082
0083     // Bizarre, pathological limits:
0084     if (variance == 0)
0085     {
0086         if (runs == expected_runs)
0087         {
0088             Real statistic = 0;
0089             Real pvalue = 1;
0090             return std::make_pair(statistic, pvalue);
0091         }
0092         else
0093         {
0094             return std::make_pair(std::numeric_limits<Real>::quiet_NaN(), Real(0));
0095         }
0096     }
0097
0098     Real sd = sqrt(variance);
0099     Real statistic = (runs - expected_runs)/sd;
0100
0101     auto normal = boost::math::normal_distribution<Real, no_promote_policy>(0,1);
0102     Real pvalue = 2*boost::math::cdf(normal, -abs(statistic));
0103     return std::make_pair(statistic, pvalue);
0104 }
0105
0106 template<class RandomAccessContainer>
0107 auto runs_above_and_below_median(RandomAccessContainer const & v)
0108 {
0109     using Real = typename RandomAccessContainer::value_type;
0110     using std::log;
0111     using std::sqrt;
0112
0113     // We have to memcpy v because the median does a partial sort,
0114     // and that would be catastrophic for the runs test.
0115     auto w = v;
0116     Real median = boost::math::statistics::median(w);
0117     return runs_above_and_below_threshold(v, median);
0118 }
0119
0120 }
0121 #endif