ceph/src/boost/boost/math/statistics/t_test.hpp

   1 //  (C) Copyright Nick Thompson 2019.
   2 //  (C) Copyright Matt Borland 2021.
   3 //  Use, modification and distribution are subject to the
   4 //  Boost Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6
   7 #ifndef BOOST_MATH_STATISTICS_T_TEST_HPP
   8 #define BOOST_MATH_STATISTICS_T_TEST_HPP
   9
  10 #include <cmath>
  11 #include <cstddef>
  12 #include <iterator>
  13 #include <utility>
  14 #include <type_traits>
  15 #include <vector>
  16 #include <stdexcept>
  17 #include <boost/math/distributions/students_t.hpp>
  18 #include <boost/math/statistics/univariate_statistics.hpp>
  19
  20 namespace boost { namespace math { namespace statistics { namespace detail {
  21
  22 template<typename ReturnType, typename T>
  23 ReturnType one_sample_t_test_impl(T sample_mean, T sample_variance, T num_samples, T assumed_mean)
  24 {
  25     using Real = typename std::tuple_element<0, ReturnType>::type;
  26     using std::sqrt;
  27     typedef boost::math::policies::policy<
  28           boost::math::policies::promote_float<false>,
  29           boost::math::policies::promote_double<false> >
  30           no_promote_policy;
  31
  32     Real test_statistic = (sample_mean - assumed_mean)/sqrt(sample_variance/num_samples);
  33     auto student = boost::math::students_t_distribution<Real, no_promote_policy>(num_samples - 1);
  34     Real pvalue;
  35     if (test_statistic > 0) {
  36         pvalue = 2*boost::math::cdf<Real>(student, -test_statistic);;
  37     }
  38     else {
  39         pvalue = 2*boost::math::cdf<Real>(student, test_statistic);
  40     }
  41     return std::make_pair(test_statistic, pvalue);
  42 }
  43
  44 template<typename ReturnType, typename ForwardIterator>
  45 ReturnType one_sample_t_test_impl(ForwardIterator begin, ForwardIterator end, typename std::iterator_traits<ForwardIterator>::value_type assumed_mean)
  46 {
  47     using Real = typename std::tuple_element<0, ReturnType>::type;
  48     std::pair<Real, Real> temp = mean_and_sample_variance(begin, end);
  49     Real mu = std::get<0>(temp);
  50     Real s_sq = std::get<1>(temp);
  51     return one_sample_t_test_impl<ReturnType>(mu, s_sq, Real(std::distance(begin, end)), Real(assumed_mean));
  52 }
  53
  54 // https://en.wikipedia.org/wiki/Student%27s_t-test#Equal_or_unequal_sample_sizes,_unequal_variances_(sX1_%3E_2sX2_or_sX2_%3E_2sX1)
  55 template<typename ReturnType, typename T>
  56 ReturnType welchs_t_test_impl(T mean_1, T variance_1, T size_1, T mean_2, T variance_2, T size_2)
  57 {
  58     using Real = typename std::tuple_element<0, ReturnType>::type;
  59     using no_promote_policy = boost::math::policies::policy<boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>>;
  60     using std::sqrt;
  61
  62     Real dof_num = (variance_1/size_1 + variance_2/size_2) * (variance_1/size_1 + variance_2/size_2);
  63     Real dof_denom = ((variance_1/size_1) * (variance_1/size_1))/(size_1 - 1) +
  64                      ((variance_2/size_2) * (variance_2/size_2))/(size_2 - 1);
  65     Real dof = dof_num / dof_denom;
  66
  67     Real s_estimator = sqrt((variance_1/size_1) + (variance_2/size_2));
  68
  69     Real test_statistic = (static_cast<Real>(mean_1) - static_cast<Real>(mean_2))/s_estimator;
  70     auto student = boost::math::students_t_distribution<Real, no_promote_policy>(dof);
  71     Real pvalue;
  72     if (test_statistic > 0)
  73     {
  74         pvalue = 2*boost::math::cdf<Real>(student, -test_statistic);;
  75     }
  76     else
  77     {
  78         pvalue = 2*boost::math::cdf<Real>(student, test_statistic);
  79     }
  80
  81     return std::make_pair(test_statistic, pvalue);
  82 }
  83
  84 // https://en.wikipedia.org/wiki/Student%27s_t-test#Equal_or_unequal_sample_sizes,_similar_variances_(1/2_%3C_sX1/sX2_%3C_2)
  85 template<typename ReturnType, typename T>
  86 ReturnType two_sample_t_test_impl(T mean_1, T variance_1, T size_1, T mean_2, T variance_2, T size_2)
  87 {
  88     using Real = typename std::tuple_element<0, ReturnType>::type;
  89     using no_promote_policy = boost::math::policies::policy<boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>>;
  90     using std::sqrt;
  91
  92     Real dof = size_1 + size_2 - 2;
  93     Real pooled_std_dev = sqrt(((size_1-1)*variance_1 + (size_2-1)*variance_2) / dof);
  94     Real test_statistic = (mean_1-mean_2) / (pooled_std_dev*sqrt(1.0/static_cast<Real>(size_1) + 1.0/static_cast<Real>(size_2)));
  95
  96     auto student = boost::math::students_t_distribution<Real, no_promote_policy>(dof);
  97     Real pvalue;
  98     if (test_statistic > 0)
  99     {
 100         pvalue = 2*boost::math::cdf<Real>(student, -test_statistic);;
 101     }
 102     else
 103     {
 104         pvalue = 2*boost::math::cdf<Real>(student, test_statistic);
 105     }
 106
 107     return std::make_pair(test_statistic, pvalue);
 108 }
 109
 110 template<typename ReturnType, typename ForwardIterator>
 111 ReturnType two_sample_t_test_impl(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2)
 112 {
 113     using Real = typename std::tuple_element<0, ReturnType>::type;
 114     using std::sqrt;
 115     auto n1 = std::distance(begin_1, end_1);
 116     auto n2 = std::distance(begin_2, end_2);
 117
 118     ReturnType temp_1 = mean_and_sample_variance(begin_1, end_1);
 119     Real mean_1 = std::get<0>(temp_1);
 120     Real variance_1 = std::get<1>(temp_1);
 121     Real std_dev_1 = sqrt(variance_1);
 122
 123     ReturnType temp_2 = mean_and_sample_variance(begin_2, end_2);
 124     Real mean_2 = std::get<0>(temp_2);
 125     Real variance_2 = std::get<1>(temp_2);
 126     Real std_dev_2 = sqrt(variance_2);
 127
 128     if(std_dev_1 > 2 * std_dev_2 || std_dev_2 > 2 * std_dev_1)
 129     {
 130         return welchs_t_test_impl<ReturnType>(mean_1, variance_1, Real(n1), mean_2, variance_2, Real(n2));
 131     }
 132     else
 133     {
 134         return two_sample_t_test_impl<ReturnType>(mean_1, variance_1, Real(n1), mean_2, variance_2, Real(n2));
 135     }
 136 }
 137
 138 // https://en.wikipedia.org/wiki/Student%27s_t-test#Dependent_t-test_for_paired_samples
 139 template<typename ReturnType, typename ForwardIterator>
 140 ReturnType paired_samples_t_test_impl(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2)
 141 {
 142     using Real = typename std::tuple_element<0, ReturnType>::type;
 143     using no_promote_policy = boost::math::policies::policy<boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>>;
 144     using std::sqrt;
 145
 146     std::vector<Real> delta;
 147     ForwardIterator it_1 = begin_1;
 148     ForwardIterator it_2 = begin_2;
 149     std::size_t n = 0;
 150     while(it_1 != end_1 && it_2 != end_2)
 151     {
 152         delta.emplace_back(static_cast<Real>(*it_1++) - static_cast<Real>(*it_2++));
 153         ++n;
 154     }
 155
 156     if(it_1 != end_1 || it_2 != end_2)
 157     {
 158         throw std::domain_error("Both sets must have the same number of values.");
 159     }
 160
 161     std::pair<Real, Real> temp = mean_and_sample_variance(delta.begin(), delta.end());
 162     Real delta_mean = std::get<0>(temp);
 163     Real delta_std_dev = sqrt(std::get<1>(temp));
 164
 165     Real test_statistic = delta_mean/(delta_std_dev/sqrt(n));
 166
 167     auto student = boost::math::students_t_distribution<Real, no_promote_policy>(n - 1);
 168     Real pvalue;
 169     if (test_statistic > 0)
 170     {
 171         pvalue = 2*boost::math::cdf<Real>(student, -test_statistic);;
 172     }
 173     else
 174     {
 175         pvalue = 2*boost::math::cdf<Real>(student, test_statistic);
 176     }
 177
 178     return std::make_pair(test_statistic, pvalue);
 179 }
 180 } // namespace detail
 181
 182 template<typename Real, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
 183 inline auto one_sample_t_test(Real sample_mean, Real sample_variance, Real num_samples, Real assumed_mean) -> std::pair<double, double>
 184 {
 185     return detail::one_sample_t_test_impl<std::pair<double, double>>(sample_mean, sample_variance, num_samples, assumed_mean);
 186 }
 187
 188 template<typename Real, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
 189 inline auto one_sample_t_test(Real sample_mean, Real sample_variance, Real num_samples, Real assumed_mean) -> std::pair<Real, Real>
 190 {
 191     return detail::one_sample_t_test_impl<std::pair<Real, Real>>(sample_mean, sample_variance, num_samples, assumed_mean);
 192 }
 193
 194 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
 195          typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
 196 inline auto one_sample_t_test(ForwardIterator begin, ForwardIterator end, Real assumed_mean) -> std::pair<double, double>
 197 {
 198     return detail::one_sample_t_test_impl<std::pair<double, double>>(begin, end, assumed_mean);
 199 }
 200
 201 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
 202          typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
 203 inline auto one_sample_t_test(ForwardIterator begin, ForwardIterator end, Real assumed_mean) -> std::pair<Real, Real>
 204 {
 205     return detail::one_sample_t_test_impl<std::pair<Real, Real>>(begin, end, assumed_mean);
 206 }
 207
 208 template<typename Container, typename Real = typename Container::value_type,
 209          typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
 210 inline auto one_sample_t_test(Container const & v, Real assumed_mean) -> std::pair<double, double>
 211 {
 212     return detail::one_sample_t_test_impl<std::pair<double, double>>(std::begin(v), std::end(v), assumed_mean);
 213 }
 214
 215 template<typename Container, typename Real = typename Container::value_type,
 216          typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
 217 inline auto one_sample_t_test(Container const & v, Real assumed_mean) -> std::pair<Real, Real>
 218 {
 219     return detail::one_sample_t_test_impl<std::pair<Real, Real>>(std::begin(v), std::end(v), assumed_mean);
 220 }
 221
 222 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
 223          typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
 224 inline auto two_sample_t_test(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2) -> std::pair<double, double>
 225 {
 226     return detail::two_sample_t_test_impl<std::pair<double, double>>(begin_1, end_1, begin_2, end_2);
 227 }
 228
 229 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
 230          typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
 231 inline auto two_sample_t_test(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2) -> std::pair<Real, Real>
 232 {
 233     return detail::two_sample_t_test_impl<std::pair<Real, Real>>(begin_1, end_1, begin_2, end_2);
 234 }
 235
 236 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
 237 inline auto two_sample_t_test(Container const & u, Container const & v) -> std::pair<double, double>
 238 {
 239     return detail::two_sample_t_test_impl<std::pair<double, double>>(std::begin(u), std::end(u), std::begin(v), std::end(v));
 240 }
 241
 242 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
 243 inline auto two_sample_t_test(Container const & u, Container const & v) -> std::pair<Real, Real>
 244 {
 245     return detail::two_sample_t_test_impl<std::pair<Real, Real>>(std::begin(u), std::end(u), std::begin(v), std::end(v));
 246 }
 247
 248 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
 249          typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
 250 inline auto paired_samples_t_test(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2) -> std::pair<double, double>
 251 {
 252     return detail::paired_samples_t_test_impl<std::pair<double, double>>(begin_1, end_1, begin_2, end_2);
 253 }
 254
 255 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
 256          typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
 257 inline auto paired_samples_t_test(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2) -> std::pair<Real, Real>
 258 {
 259     return detail::paired_samples_t_test_impl<std::pair<Real, Real>>(begin_1, end_1, begin_2, end_2);
 260 }
 261
 262 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
 263 inline auto paired_samples_t_test(Container const & u, Container const & v) -> std::pair<double, double>
 264 {
 265     return detail::paired_samples_t_test_impl<std::pair<double, double>>(std::begin(u), std::end(u), std::begin(v), std::end(v));
 266 }
 267
 268 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
 269 inline auto paired_samples_t_test(Container const & u, Container const & v) -> std::pair<Real, Real>
 270 {
 271     return detail::paired_samples_t_test_impl<std::pair<Real, Real>>(std::begin(u), std::end(u), std::begin(v), std::end(v));
 272 }
 273
 274 }}} // namespace boost::math::statistics
 275 #endif