276 lines · plain
1// (C) Copyright Nick Thompson 2019.2// (C) Copyright Matt Borland 2021.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_STATISTICS_T_TEST_HPP8#define BOOST_MATH_STATISTICS_T_TEST_HPP9 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 20namespace boost { namespace math { namespace statistics { namespace detail {21 22template<typename ReturnType, typename T>23ReturnType 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 44template<typename ReturnType, typename ForwardIterator>45ReturnType 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)55template<typename ReturnType, typename T>56ReturnType 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)85template<typename ReturnType, typename T>86ReturnType 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 110template<typename ReturnType, typename ForwardIterator>111ReturnType 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 else133 {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_samples139template<typename ReturnType, typename ForwardIterator>140ReturnType 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 detail181 182template<typename Real, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>183inline 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 188template<typename Real, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>189inline 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 194template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, 195 typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>196inline 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 201template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, 202 typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>203inline 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 208template<typename Container, typename Real = typename Container::value_type,209 typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>210inline 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 215template<typename Container, typename Real = typename Container::value_type,216 typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>217inline 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 222template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, 223 typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>224inline 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 229template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, 230 typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>231inline 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 236template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>237inline 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 242template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>243inline 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 248template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, 249 typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>250inline 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 255template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, 256 typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>257inline 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 262template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>263inline 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 268template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>269inline 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::statistics275#endif276