From 2adaa0334b60db81bf46e7bf47f1d39eeeb2f346 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Sun, 13 Oct 2019 10:02:02 -0400 Subject: [PATCH] Runs test [CI SKIP] --- include/boost/math/statistics/runs_test.hpp | 64 +++++++++++++++++++++ 1 file changed, 64 insertions(+) create mode 100644 include/boost/math/statistics/runs_test.hpp diff --git a/include/boost/math/statistics/runs_test.hpp b/include/boost/math/statistics/runs_test.hpp new file mode 100644 index 000000000..ea32a7f1d --- /dev/null +++ b/include/boost/math/statistics/runs_test.hpp @@ -0,0 +1,64 @@ +/* + * Copyright Nick Thompson, 2019 + * Use, modification and distribution are subject to the + * Boost Software License, Version 1.0. (See accompanying file + * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + */ + +#ifndef BOOST_MATH_STATISTICS_RUNS_TEST_HPP +#define BOOST_MATH_STATISTICS_RUNS_TEST_HPP + +#include +#include +#include + +namespace boost::math::statistics { + +template +auto runs_above_threshold(RandomAccessContainer const & v, + typename RandomAccessContainer::value_type threshold) +{ + typedef boost::math::policies::policy< + boost::math::policies::promote_float, + boost::math::policies::promote_double > + no_promote_policy; + + decltype(v.size()) nabove = 0; + decltype(v.size()) nbelow = 0; + + for (auto const & x : v) { + if (x > threshold) { + ++nabove; + } + if (x < threshold) { + ++nbelow; + } + } + Real n = nabove + nbelow; + + Real mu = Real(1) + Real(2*nabove*nbelow)/Real(n); + Real var = 2*nabove*nbelow*(2*nabove*nbelow-n)/Real(n*n*(n-1)); +} + +template +auto runs_above_median(RandomAccessContainer const & v) +{ + using Real = typename RandomAccessContainer::value_type; + using std::log; + using std::sqrt; + + Real median; + + { + // We have to memcpy v because the median does a partial sort, + // and that would be catastrophic for the runs test. + auto w = v; + median = boost::math::statistics::median(w); + } + + + return runs_above_threshold(v, median); +} + +} +#endif