math/statistics/t_test.hpp

0001 //  (C) Copyright Nick Thompson 2019.
0002 //  (C) Copyright Matt Borland 2021.
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 #ifndef BOOST_MATH_STATISTICS_T_TEST_HPP
0008 #define BOOST_MATH_STATISTICS_T_TEST_HPP
0009
0010 #include <cmath>
0011 #include <cstddef>
0012 #include <iterator>
0013 #include <utility>
0014 #include <type_traits>
0015 #include <vector>
0016 #include <stdexcept>
0017 #include <boost/math/distributions/students_t.hpp>
0018 #include <boost/math/statistics/univariate_statistics.hpp>
0019
0020 namespace boost { namespace math { namespace statistics { namespace detail {
0021
0022 template<typename ReturnType, typename T>
0023 ReturnType one_sample_t_test_impl(T sample_mean, T sample_variance, T num_samples, T assumed_mean)
0024 {
0025     using Real = typename std::tuple_element<0, ReturnType>::type;
0026     using std::sqrt;
0027     typedef boost::math::policies::policy<
0028           boost::math::policies::promote_float<false>,
0029           boost::math::policies::promote_double<false> >
0030           no_promote_policy;
0031
0032     Real test_statistic = (sample_mean - assumed_mean)/sqrt(sample_variance/num_samples);
0033     auto student = boost::math::students_t_distribution<Real, no_promote_policy>(num_samples - 1);
0034     Real pvalue;
0035     if (test_statistic > 0) {
0036         pvalue = 2*boost::math::cdf<Real>(student, -test_statistic);;
0037     }
0038     else {
0039         pvalue = 2*boost::math::cdf<Real>(student, test_statistic);
0040     }
0041     return std::make_pair(test_statistic, pvalue);
0042 }
0043
0044 template<typename ReturnType, typename ForwardIterator>
0045 ReturnType one_sample_t_test_impl(ForwardIterator begin, ForwardIterator end, typename std::iterator_traits<ForwardIterator>::value_type assumed_mean)
0046 {
0047     using Real = typename std::tuple_element<0, ReturnType>::type;
0048     std::pair<Real, Real> temp = mean_and_sample_variance(begin, end);
0049     Real mu = std::get<0>(temp);
0050     Real s_sq = std::get<1>(temp);
0051     return one_sample_t_test_impl<ReturnType>(mu, s_sq, Real(std::distance(begin, end)), Real(assumed_mean));
0052 }
0053
0054 // https://en.wikipedia.org/wiki/Student%27s_t-test#Equal_or_unequal_sample_sizes,_unequal_variances_(sX1_%3E_2sX2_or_sX2_%3E_2sX1)
0055 template<typename ReturnType, typename T>
0056 ReturnType welchs_t_test_impl(T mean_1, T variance_1, T size_1, T mean_2, T variance_2, T size_2)
0057 {
0058     using Real = typename std::tuple_element<0, ReturnType>::type;
0059     using no_promote_policy = boost::math::policies::policy<boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>>;
0060     using std::sqrt;
0061
0062     Real dof_num = (variance_1/size_1 + variance_2/size_2) * (variance_1/size_1 + variance_2/size_2);
0063     Real dof_denom = ((variance_1/size_1) * (variance_1/size_1))/(size_1 - 1) +
0064                      ((variance_2/size_2) * (variance_2/size_2))/(size_2 - 1);
0065     Real dof = dof_num / dof_denom;
0066
0067     Real s_estimator = sqrt((variance_1/size_1) + (variance_2/size_2));
0068
0069     Real test_statistic = (static_cast<Real>(mean_1) - static_cast<Real>(mean_2))/s_estimator;
0070     auto student = boost::math::students_t_distribution<Real, no_promote_policy>(dof);
0071     Real pvalue;
0072     if (test_statistic > 0)
0073     {
0074         pvalue = 2*boost::math::cdf<Real>(student, -test_statistic);;
0075     }
0076     else
0077     {
0078         pvalue = 2*boost::math::cdf<Real>(student, test_statistic);
0079     }
0080
0081     return std::make_pair(test_statistic, pvalue);
0082 }
0083
0084 // https://en.wikipedia.org/wiki/Student%27s_t-test#Equal_or_unequal_sample_sizes,_similar_variances_(1/2_%3C_sX1/sX2_%3C_2)
0085 template<typename ReturnType, typename T>
0086 ReturnType two_sample_t_test_impl(T mean_1, T variance_1, T size_1, T mean_2, T variance_2, T size_2)
0087 {
0088     using Real = typename std::tuple_element<0, ReturnType>::type;
0089     using no_promote_policy = boost::math::policies::policy<boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>>;
0090     using std::sqrt;
0091
0092     Real dof = size_1 + size_2 - 2;
0093     Real pooled_std_dev = sqrt(((size_1-1)*variance_1 + (size_2-1)*variance_2) / dof);
0094     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)));
0095
0096     auto student = boost::math::students_t_distribution<Real, no_promote_policy>(dof);
0097     Real pvalue;
0098     if (test_statistic > 0)
0099     {
0100         pvalue = 2*boost::math::cdf<Real>(student, -test_statistic);;
0101     }
0102     else
0103     {
0104         pvalue = 2*boost::math::cdf<Real>(student, test_statistic);
0105     }
0106
0107     return std::make_pair(test_statistic, pvalue);
0108 }
0109
0110 template<typename ReturnType, typename ForwardIterator>
0111 ReturnType two_sample_t_test_impl(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2)
0112 {
0113     using Real = typename std::tuple_element<0, ReturnType>::type;
0114     using std::sqrt;
0115     auto n1 = std::distance(begin_1, end_1);
0116     auto n2 = std::distance(begin_2, end_2);
0117
0118     ReturnType temp_1 = mean_and_sample_variance(begin_1, end_1);
0119     Real mean_1 = std::get<0>(temp_1);
0120     Real variance_1 = std::get<1>(temp_1);
0121     Real std_dev_1 = sqrt(variance_1);
0122
0123     ReturnType temp_2 = mean_and_sample_variance(begin_2, end_2);
0124     Real mean_2 = std::get<0>(temp_2);
0125     Real variance_2 = std::get<1>(temp_2);
0126     Real std_dev_2 = sqrt(variance_2);
0127
0128     if(std_dev_1 > 2 * std_dev_2 || std_dev_2 > 2 * std_dev_1)
0129     {
0130         return welchs_t_test_impl<ReturnType>(mean_1, variance_1, Real(n1), mean_2, variance_2, Real(n2));
0131     }
0132     else
0133     {
0134         return two_sample_t_test_impl<ReturnType>(mean_1, variance_1, Real(n1), mean_2, variance_2, Real(n2));
0135     }
0136 }
0137
0138 // https://en.wikipedia.org/wiki/Student%27s_t-test#Dependent_t-test_for_paired_samples
0139 template<typename ReturnType, typename ForwardIterator>
0140 ReturnType paired_samples_t_test_impl(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2)
0141 {
0142     using Real = typename std::tuple_element<0, ReturnType>::type;
0143     using no_promote_policy = boost::math::policies::policy<boost::math::policies::promote_float<false>, boost::math::policies::promote_double<false>>;
0144     using std::sqrt;
0145
0146     std::vector<Real> delta;
0147     ForwardIterator it_1 = begin_1;
0148     ForwardIterator it_2 = begin_2;
0149     std::size_t n = 0;
0150     while(it_1 != end_1 && it_2 != end_2)
0151     {
0152         delta.emplace_back(static_cast<Real>(*it_1++) - static_cast<Real>(*it_2++));
0153         ++n;
0154     }
0155
0156     if(it_1 != end_1 || it_2 != end_2)
0157     {
0158         throw std::domain_error("Both sets must have the same number of values.");
0159     }
0160
0161     std::pair<Real, Real> temp = mean_and_sample_variance(delta.begin(), delta.end());
0162     Real delta_mean = std::get<0>(temp);
0163     Real delta_std_dev = sqrt(std::get<1>(temp));
0164
0165     Real test_statistic = delta_mean/(delta_std_dev/sqrt(n));
0166
0167     auto student = boost::math::students_t_distribution<Real, no_promote_policy>(n - 1);
0168     Real pvalue;
0169     if (test_statistic > 0)
0170     {
0171         pvalue = 2*boost::math::cdf<Real>(student, -test_statistic);;
0172     }
0173     else
0174     {
0175         pvalue = 2*boost::math::cdf<Real>(student, test_statistic);
0176     }
0177
0178     return std::make_pair(test_statistic, pvalue);
0179 }
0180 } // namespace detail
0181
0182 template<typename Real, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
0183 inline auto one_sample_t_test(Real sample_mean, Real sample_variance, Real num_samples, Real assumed_mean) -> std::pair<double, double>
0184 {
0185     return detail::one_sample_t_test_impl<std::pair<double, double>>(sample_mean, sample_variance, num_samples, assumed_mean);
0186 }
0187
0188 template<typename Real, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
0189 inline auto one_sample_t_test(Real sample_mean, Real sample_variance, Real num_samples, Real assumed_mean) -> std::pair<Real, Real>
0190 {
0191     return detail::one_sample_t_test_impl<std::pair<Real, Real>>(sample_mean, sample_variance, num_samples, assumed_mean);
0192 }
0193
0194 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
0195          typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
0196 inline auto one_sample_t_test(ForwardIterator begin, ForwardIterator end, Real assumed_mean) -> std::pair<double, double>
0197 {
0198     return detail::one_sample_t_test_impl<std::pair<double, double>>(begin, end, assumed_mean);
0199 }
0200
0201 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
0202          typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
0203 inline auto one_sample_t_test(ForwardIterator begin, ForwardIterator end, Real assumed_mean) -> std::pair<Real, Real>
0204 {
0205     return detail::one_sample_t_test_impl<std::pair<Real, Real>>(begin, end, assumed_mean);
0206 }
0207
0208 template<typename Container, typename Real = typename Container::value_type,
0209          typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
0210 inline auto one_sample_t_test(Container const & v, Real assumed_mean) -> std::pair<double, double>
0211 {
0212     return detail::one_sample_t_test_impl<std::pair<double, double>>(std::begin(v), std::end(v), assumed_mean);
0213 }
0214
0215 template<typename Container, typename Real = typename Container::value_type,
0216          typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
0217 inline auto one_sample_t_test(Container const & v, Real assumed_mean) -> std::pair<Real, Real>
0218 {
0219     return detail::one_sample_t_test_impl<std::pair<Real, Real>>(std::begin(v), std::end(v), assumed_mean);
0220 }
0221
0222 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
0223          typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
0224 inline auto two_sample_t_test(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2) -> std::pair<double, double>
0225 {
0226     return detail::two_sample_t_test_impl<std::pair<double, double>>(begin_1, end_1, begin_2, end_2);
0227 }
0228
0229 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
0230          typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
0231 inline auto two_sample_t_test(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2) -> std::pair<Real, Real>
0232 {
0233     return detail::two_sample_t_test_impl<std::pair<Real, Real>>(begin_1, end_1, begin_2, end_2);
0234 }
0235
0236 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
0237 inline auto two_sample_t_test(Container const & u, Container const & v) -> std::pair<double, double>
0238 {
0239     return detail::two_sample_t_test_impl<std::pair<double, double>>(std::begin(u), std::end(u), std::begin(v), std::end(v));
0240 }
0241
0242 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
0243 inline auto two_sample_t_test(Container const & u, Container const & v) -> std::pair<Real, Real>
0244 {
0245     return detail::two_sample_t_test_impl<std::pair<Real, Real>>(std::begin(u), std::end(u), std::begin(v), std::end(v));
0246 }
0247
0248 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
0249          typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
0250 inline auto paired_samples_t_test(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2) -> std::pair<double, double>
0251 {
0252     return detail::paired_samples_t_test_impl<std::pair<double, double>>(begin_1, end_1, begin_2, end_2);
0253 }
0254
0255 template<typename ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type,
0256          typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
0257 inline auto paired_samples_t_test(ForwardIterator begin_1, ForwardIterator end_1, ForwardIterator begin_2, ForwardIterator end_2) -> std::pair<Real, Real>
0258 {
0259     return detail::paired_samples_t_test_impl<std::pair<Real, Real>>(begin_1, end_1, begin_2, end_2);
0260 }
0261
0262 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<std::is_integral<Real>::value, bool>::type = true>
0263 inline auto paired_samples_t_test(Container const & u, Container const & v) -> std::pair<double, double>
0264 {
0265     return detail::paired_samples_t_test_impl<std::pair<double, double>>(std::begin(u), std::end(u), std::begin(v), std::end(v));
0266 }
0267
0268 template<typename Container, typename Real = typename Container::value_type, typename std::enable_if<!std::is_integral<Real>::value, bool>::type = true>
0269 inline auto paired_samples_t_test(Container const & u, Container const & v) -> std::pair<Real, Real>
0270 {
0271     return detail::paired_samples_t_test_impl<std::pair<Real, Real>>(std::begin(u), std::end(u), std::begin(v), std::end(v));
0272 }
0273
0274 }}} // namespace boost::math::statistics
0275 #endif