t_test.hpp 13 KB

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