]>
Commit | Line | Data |
---|---|---|
92f5a8d4 | 1 | // (C) Copyright Nick Thompson 2018. |
1e59de90 | 2 | // (C) Copyright Matt Borland 2020. |
92f5a8d4 TL |
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_UNIVARIATE_STATISTICS_HPP | |
8 | #define BOOST_MATH_STATISTICS_UNIVARIATE_STATISTICS_HPP | |
9 | ||
1e59de90 TL |
10 | #include <boost/math/statistics/detail/single_pass.hpp> |
11 | #include <boost/math/tools/config.hpp> | |
12 | #include <boost/math/tools/assert.hpp> | |
92f5a8d4 TL |
13 | #include <algorithm> |
14 | #include <iterator> | |
f67539c2 TL |
15 | #include <tuple> |
16 | #include <cmath> | |
20effc67 | 17 | #include <vector> |
1e59de90 TL |
18 | #include <type_traits> |
19 | #include <utility> | |
20 | #include <numeric> | |
21 | #include <list> | |
22 | ||
23 | // Support compilers with P0024R2 implemented without linking TBB | |
24 | // https://en.cppreference.com/w/cpp/compiler_support | |
25 | #if !defined(BOOST_NO_CXX17_HDR_EXECUTION) && defined(BOOST_HAS_THREADS) | |
26 | #include <execution> | |
92f5a8d4 TL |
27 | |
28 | namespace boost::math::statistics { | |
29 | ||
1e59de90 TL |
30 | template<class ExecutionPolicy, class ForwardIterator> |
31 | inline auto mean(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last) | |
92f5a8d4 TL |
32 | { |
33 | using Real = typename std::iterator_traits<ForwardIterator>::value_type; | |
1e59de90 TL |
34 | BOOST_MATH_ASSERT_MSG(first != last, "At least one sample is required to compute the mean."); |
35 | ||
36 | if constexpr (std::is_integral_v<Real>) | |
37 | { | |
38 | if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>) | |
39 | { | |
40 | return detail::mean_sequential_impl<double>(first, last); | |
92f5a8d4 | 41 | } |
1e59de90 | 42 | else |
92f5a8d4 | 43 | { |
1e59de90 | 44 | return std::reduce(exec, first, last, 0.0) / std::distance(first, last); |
92f5a8d4 | 45 | } |
92f5a8d4 TL |
46 | } |
47 | else | |
48 | { | |
1e59de90 | 49 | if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>) |
92f5a8d4 | 50 | { |
1e59de90 TL |
51 | return detail::mean_sequential_impl<Real>(first, last); |
52 | } | |
53 | else | |
54 | { | |
55 | return std::reduce(exec, first, last, Real(0.0)) / Real(std::distance(first, last)); | |
92f5a8d4 | 56 | } |
92f5a8d4 TL |
57 | } |
58 | } | |
59 | ||
1e59de90 TL |
60 | template<class ExecutionPolicy, class Container> |
61 | inline auto mean(ExecutionPolicy&& exec, Container const & v) | |
62 | { | |
63 | return mean(exec, std::cbegin(v), std::cend(v)); | |
64 | } | |
65 | ||
66 | template<class ForwardIterator> | |
67 | inline auto mean(ForwardIterator first, ForwardIterator last) | |
68 | { | |
69 | return mean(std::execution::seq, first, last); | |
70 | } | |
71 | ||
92f5a8d4 TL |
72 | template<class Container> |
73 | inline auto mean(Container const & v) | |
74 | { | |
1e59de90 | 75 | return mean(std::execution::seq, std::cbegin(v), std::cend(v)); |
92f5a8d4 TL |
76 | } |
77 | ||
1e59de90 TL |
78 | template<class ExecutionPolicy, class ForwardIterator> |
79 | inline auto variance(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last) | |
92f5a8d4 TL |
80 | { |
81 | using Real = typename std::iterator_traits<ForwardIterator>::value_type; | |
1e59de90 TL |
82 | |
83 | if constexpr (std::is_integral_v<Real>) | |
84 | { | |
85 | if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>) | |
92f5a8d4 | 86 | { |
1e59de90 TL |
87 | return std::get<2>(detail::variance_sequential_impl<std::tuple<double, double, double, double>>(first, last)); |
88 | } | |
89 | else | |
90 | { | |
91 | const auto results = detail::first_four_moments_parallel_impl<std::tuple<double, double, double, double, double>>(first, last); | |
92 | return std::get<1>(results) / std::get<4>(results); | |
92f5a8d4 | 93 | } |
92f5a8d4 TL |
94 | } |
95 | else | |
96 | { | |
1e59de90 TL |
97 | if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>) |
98 | { | |
99 | return std::get<2>(detail::variance_sequential_impl<std::tuple<Real, Real, Real, Real>>(first, last)); | |
100 | } | |
101 | else | |
92f5a8d4 | 102 | { |
1e59de90 TL |
103 | const auto results = detail::first_four_moments_parallel_impl<std::tuple<Real, Real, Real, Real, Real>>(first, last); |
104 | return std::get<1>(results) / std::get<4>(results); | |
92f5a8d4 | 105 | } |
92f5a8d4 TL |
106 | } |
107 | } | |
108 | ||
1e59de90 TL |
109 | template<class ExecutionPolicy, class Container> |
110 | inline auto variance(ExecutionPolicy&& exec, Container const & v) | |
111 | { | |
112 | return variance(exec, std::cbegin(v), std::cend(v)); | |
113 | } | |
114 | ||
115 | template<class ForwardIterator> | |
116 | inline auto variance(ForwardIterator first, ForwardIterator last) | |
117 | { | |
118 | return variance(std::execution::seq, first, last); | |
119 | } | |
120 | ||
92f5a8d4 TL |
121 | template<class Container> |
122 | inline auto variance(Container const & v) | |
123 | { | |
1e59de90 TL |
124 | return variance(std::execution::seq, std::cbegin(v), std::cend(v)); |
125 | } | |
126 | ||
127 | template<class ExecutionPolicy, class ForwardIterator> | |
128 | inline auto sample_variance(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last) | |
129 | { | |
130 | const auto n = std::distance(first, last); | |
131 | BOOST_MATH_ASSERT_MSG(n > 1, "At least two samples are required to compute the sample variance."); | |
132 | return n*variance(exec, first, last)/(n-1); | |
133 | } | |
134 | ||
135 | template<class ExecutionPolicy, class Container> | |
136 | inline auto sample_variance(ExecutionPolicy&& exec, Container const & v) | |
137 | { | |
138 | return sample_variance(exec, std::cbegin(v), std::cend(v)); | |
92f5a8d4 TL |
139 | } |
140 | ||
141 | template<class ForwardIterator> | |
1e59de90 | 142 | inline auto sample_variance(ForwardIterator first, ForwardIterator last) |
92f5a8d4 | 143 | { |
1e59de90 | 144 | return sample_variance(std::execution::seq, first, last); |
92f5a8d4 TL |
145 | } |
146 | ||
147 | template<class Container> | |
148 | inline auto sample_variance(Container const & v) | |
149 | { | |
1e59de90 | 150 | return sample_variance(std::execution::seq, std::cbegin(v), std::cend(v)); |
92f5a8d4 TL |
151 | } |
152 | ||
1e59de90 TL |
153 | template<class ExecutionPolicy, class ForwardIterator> |
154 | inline auto mean_and_sample_variance(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last) | |
92f5a8d4 TL |
155 | { |
156 | using Real = typename std::iterator_traits<ForwardIterator>::value_type; | |
1e59de90 TL |
157 | |
158 | if constexpr (std::is_integral_v<Real>) | |
159 | { | |
160 | if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>) | |
92f5a8d4 | 161 | { |
1e59de90 TL |
162 | const auto results = detail::variance_sequential_impl<std::tuple<double, double, double, double>>(first, last); |
163 | return std::make_pair(std::get<0>(results), std::get<2>(results)*std::get<3>(results)/(std::get<3>(results)-1.0)); | |
164 | } | |
165 | else | |
166 | { | |
167 | const auto results = detail::first_four_moments_parallel_impl<std::tuple<double, double, double, double, double>>(first, last); | |
168 | return std::make_pair(std::get<0>(results), std::get<1>(results) / (std::get<4>(results)-1.0)); | |
92f5a8d4 | 169 | } |
92f5a8d4 TL |
170 | } |
171 | else | |
172 | { | |
1e59de90 TL |
173 | if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>) |
174 | { | |
175 | const auto results = detail::variance_sequential_impl<std::tuple<Real, Real, Real, Real>>(first, last); | |
176 | return std::make_pair(std::get<0>(results), std::get<2>(results)*std::get<3>(results)/(std::get<3>(results)-Real(1))); | |
177 | } | |
178 | else | |
92f5a8d4 | 179 | { |
1e59de90 TL |
180 | const auto results = detail::first_four_moments_parallel_impl<std::tuple<Real, Real, Real, Real, Real>>(first, last); |
181 | return std::make_pair(std::get<0>(results), std::get<1>(results) / (std::get<4>(results)-Real(1))); | |
92f5a8d4 | 182 | } |
92f5a8d4 TL |
183 | } |
184 | } | |
185 | ||
1e59de90 TL |
186 | template<class ExecutionPolicy, class Container> |
187 | inline auto mean_and_sample_variance(ExecutionPolicy&& exec, Container const & v) | |
92f5a8d4 | 188 | { |
1e59de90 | 189 | return mean_and_sample_variance(exec, std::cbegin(v), std::cend(v)); |
92f5a8d4 TL |
190 | } |
191 | ||
92f5a8d4 | 192 | template<class ForwardIterator> |
1e59de90 TL |
193 | inline auto mean_and_sample_variance(ForwardIterator first, ForwardIterator last) |
194 | { | |
195 | return mean_and_sample_variance(std::execution::seq, first, last); | |
196 | } | |
197 | ||
198 | template<class Container> | |
199 | inline auto mean_and_sample_variance(Container const & v) | |
200 | { | |
201 | return mean_and_sample_variance(std::execution::seq, std::cbegin(v), std::cend(v)); | |
202 | } | |
203 | ||
204 | template<class ExecutionPolicy, class ForwardIterator> | |
205 | inline auto first_four_moments(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last) | |
92f5a8d4 TL |
206 | { |
207 | using Real = typename std::iterator_traits<ForwardIterator>::value_type; | |
1e59de90 TL |
208 | |
209 | if constexpr (std::is_integral_v<Real>) | |
210 | { | |
211 | if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>) | |
92f5a8d4 | 212 | { |
1e59de90 TL |
213 | const auto results = detail::first_four_moments_sequential_impl<std::tuple<double, double, double, double, double>>(first, last); |
214 | return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results), | |
215 | std::get<3>(results) / std::get<4>(results)); | |
92f5a8d4 | 216 | } |
1e59de90 | 217 | else |
92f5a8d4 | 218 | { |
1e59de90 TL |
219 | const auto results = detail::first_four_moments_parallel_impl<std::tuple<double, double, double, double, double>>(first, last); |
220 | return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results), | |
221 | std::get<3>(results) / std::get<4>(results)); | |
92f5a8d4 | 222 | } |
92f5a8d4 TL |
223 | } |
224 | else | |
225 | { | |
1e59de90 | 226 | if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>) |
92f5a8d4 | 227 | { |
1e59de90 TL |
228 | const auto results = detail::first_four_moments_sequential_impl<std::tuple<Real, Real, Real, Real, Real>>(first, last); |
229 | return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results), | |
230 | std::get<3>(results) / std::get<4>(results)); | |
92f5a8d4 | 231 | } |
1e59de90 | 232 | else |
92f5a8d4 | 233 | { |
1e59de90 TL |
234 | const auto results = detail::first_four_moments_parallel_impl<std::tuple<Real, Real, Real, Real, Real>>(first, last); |
235 | return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results), | |
236 | std::get<3>(results) / std::get<4>(results)); | |
92f5a8d4 | 237 | } |
92f5a8d4 TL |
238 | } |
239 | } | |
240 | ||
1e59de90 TL |
241 | template<class ExecutionPolicy, class Container> |
242 | inline auto first_four_moments(ExecutionPolicy&& exec, Container const & v) | |
243 | { | |
244 | return first_four_moments(exec, std::cbegin(v), std::cend(v)); | |
245 | } | |
246 | ||
247 | template<class ForwardIterator> | |
248 | inline auto first_four_moments(ForwardIterator first, ForwardIterator last) | |
249 | { | |
250 | return first_four_moments(std::execution::seq, first, last); | |
251 | } | |
252 | ||
92f5a8d4 | 253 | template<class Container> |
1e59de90 | 254 | inline auto first_four_moments(Container const & v) |
92f5a8d4 | 255 | { |
1e59de90 | 256 | return first_four_moments(std::execution::seq, std::cbegin(v), std::cend(v)); |
92f5a8d4 TL |
257 | } |
258 | ||
92f5a8d4 | 259 | // https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf |
1e59de90 TL |
260 | template<class ExecutionPolicy, class ForwardIterator> |
261 | inline auto skewness(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last) | |
92f5a8d4 TL |
262 | { |
263 | using Real = typename std::iterator_traits<ForwardIterator>::value_type; | |
1e59de90 TL |
264 | using std::sqrt; |
265 | ||
266 | if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>) | |
267 | { | |
268 | if constexpr (std::is_integral_v<Real>) | |
92f5a8d4 | 269 | { |
1e59de90 TL |
270 | return detail::skewness_sequential_impl<double>(first, last); |
271 | } | |
272 | else | |
273 | { | |
274 | return detail::skewness_sequential_impl<Real>(first, last); | |
92f5a8d4 | 275 | } |
92f5a8d4 | 276 | } |
1e59de90 | 277 | else |
92f5a8d4 | 278 | { |
1e59de90 TL |
279 | const auto [M1, M2, M3, M4] = first_four_moments(exec, first, last); |
280 | const auto n = std::distance(first, last); | |
281 | const auto var = M2/(n-1); | |
282 | ||
283 | if (M2 == 0) | |
92f5a8d4 | 284 | { |
1e59de90 TL |
285 | // The limit is technically undefined, but the interpretation here is clear: |
286 | // A constant dataset has no skewness. | |
287 | if constexpr (std::is_integral_v<Real>) | |
288 | { | |
289 | return double(0); | |
290 | } | |
291 | else | |
292 | { | |
293 | return Real(0); | |
294 | } | |
295 | } | |
296 | else | |
297 | { | |
298 | return M3/(M2*sqrt(var)) / Real(2); | |
92f5a8d4 | 299 | } |
92f5a8d4 TL |
300 | } |
301 | } | |
302 | ||
1e59de90 TL |
303 | template<class ExecutionPolicy, class Container> |
304 | inline auto skewness(ExecutionPolicy&& exec, Container & v) | |
92f5a8d4 | 305 | { |
1e59de90 | 306 | return skewness(exec, std::cbegin(v), std::cend(v)); |
92f5a8d4 TL |
307 | } |
308 | ||
1e59de90 TL |
309 | template<class ForwardIterator> |
310 | inline auto skewness(ForwardIterator first, ForwardIterator last) | |
311 | { | |
312 | return skewness(std::execution::seq, first, last); | |
313 | } | |
314 | ||
315 | template<class Container> | |
316 | inline auto skewness(Container const & v) | |
317 | { | |
318 | return skewness(std::execution::seq, std::cbegin(v), std::cend(v)); | |
319 | } | |
92f5a8d4 TL |
320 | |
321 | // Follows equation 1.6 of: | |
322 | // https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf | |
1e59de90 TL |
323 | template<class ExecutionPolicy, class ForwardIterator> |
324 | inline auto kurtosis(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last) | |
92f5a8d4 | 325 | { |
1e59de90 | 326 | const auto [M1, M2, M3, M4] = first_four_moments(exec, first, last); |
92f5a8d4 TL |
327 | if (M2 == 0) |
328 | { | |
329 | return M2; | |
330 | } | |
331 | return M4/(M2*M2); | |
332 | } | |
333 | ||
1e59de90 TL |
334 | template<class ExecutionPolicy, class Container> |
335 | inline auto kurtosis(ExecutionPolicy&& exec, Container const & v) | |
336 | { | |
337 | return kurtosis(exec, std::cbegin(v), std::cend(v)); | |
338 | } | |
339 | ||
340 | template<class ForwardIterator> | |
341 | inline auto kurtosis(ForwardIterator first, ForwardIterator last) | |
342 | { | |
343 | return kurtosis(std::execution::seq, first, last); | |
344 | } | |
345 | ||
92f5a8d4 TL |
346 | template<class Container> |
347 | inline auto kurtosis(Container const & v) | |
348 | { | |
1e59de90 TL |
349 | return kurtosis(std::execution::seq, std::cbegin(v), std::cend(v)); |
350 | } | |
351 | ||
352 | template<class ExecutionPolicy, class ForwardIterator> | |
353 | inline auto excess_kurtosis(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last) | |
354 | { | |
355 | return kurtosis(exec, first, last) - 3; | |
356 | } | |
357 | ||
358 | template<class ExecutionPolicy, class Container> | |
359 | inline auto excess_kurtosis(ExecutionPolicy&& exec, Container const & v) | |
360 | { | |
361 | return excess_kurtosis(exec, std::cbegin(v), std::cend(v)); | |
92f5a8d4 TL |
362 | } |
363 | ||
364 | template<class ForwardIterator> | |
1e59de90 | 365 | inline auto excess_kurtosis(ForwardIterator first, ForwardIterator last) |
92f5a8d4 | 366 | { |
1e59de90 | 367 | return excess_kurtosis(std::execution::seq, first, last); |
92f5a8d4 TL |
368 | } |
369 | ||
370 | template<class Container> | |
371 | inline auto excess_kurtosis(Container const & v) | |
372 | { | |
1e59de90 | 373 | return excess_kurtosis(std::execution::seq, std::cbegin(v), std::cend(v)); |
92f5a8d4 TL |
374 | } |
375 | ||
376 | ||
1e59de90 TL |
377 | template<class ExecutionPolicy, class RandomAccessIterator> |
378 | auto median(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last) | |
92f5a8d4 | 379 | { |
1e59de90 TL |
380 | const auto num_elems = std::distance(first, last); |
381 | BOOST_MATH_ASSERT_MSG(num_elems > 0, "The median of a zero length vector is undefined."); | |
92f5a8d4 TL |
382 | if (num_elems & 1) |
383 | { | |
384 | auto middle = first + (num_elems - 1)/2; | |
1e59de90 | 385 | std::nth_element(exec, first, middle, last); |
92f5a8d4 TL |
386 | return *middle; |
387 | } | |
388 | else | |
389 | { | |
390 | auto middle = first + num_elems/2 - 1; | |
1e59de90 TL |
391 | std::nth_element(exec, first, middle, last); |
392 | std::nth_element(exec, middle, middle+1, last); | |
92f5a8d4 TL |
393 | return (*middle + *(middle+1))/2; |
394 | } | |
395 | } | |
396 | ||
397 | ||
1e59de90 TL |
398 | template<class ExecutionPolicy, class RandomAccessContainer> |
399 | inline auto median(ExecutionPolicy&& exec, RandomAccessContainer & v) | |
400 | { | |
401 | return median(exec, std::begin(v), std::end(v)); | |
402 | } | |
403 | ||
404 | template<class RandomAccessIterator> | |
405 | inline auto median(RandomAccessIterator first, RandomAccessIterator last) | |
406 | { | |
407 | return median(std::execution::seq, first, last); | |
408 | } | |
409 | ||
92f5a8d4 TL |
410 | template<class RandomAccessContainer> |
411 | inline auto median(RandomAccessContainer & v) | |
412 | { | |
1e59de90 | 413 | return median(std::execution::seq, std::begin(v), std::end(v)); |
92f5a8d4 TL |
414 | } |
415 | ||
1e59de90 TL |
416 | #if 0 |
417 | // | |
418 | // Parallel gini calculation is curently broken, see: | |
419 | // https://github.com/boostorg/math/issues/585 | |
420 | // We will fix this at a later date, for now just use a serial implementation: | |
421 | // | |
422 | template<class ExecutionPolicy, class RandomAccessIterator> | |
423 | inline auto gini_coefficient(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last) | |
92f5a8d4 TL |
424 | { |
425 | using Real = typename std::iterator_traits<RandomAccessIterator>::value_type; | |
92f5a8d4 | 426 | |
1e59de90 | 427 | if(!std::is_sorted(exec, first, last)) |
92f5a8d4 | 428 | { |
1e59de90 TL |
429 | std::sort(exec, first, last); |
430 | } | |
92f5a8d4 | 431 | |
1e59de90 TL |
432 | if constexpr (std::is_same_v<std::remove_reference_t<decltype(exec)>, decltype(std::execution::seq)>) |
433 | { | |
434 | if constexpr (std::is_integral_v<Real>) | |
92f5a8d4 | 435 | { |
1e59de90 | 436 | return detail::gini_coefficient_sequential_impl<double>(first, last); |
92f5a8d4 | 437 | } |
1e59de90 TL |
438 | else |
439 | { | |
440 | return detail::gini_coefficient_sequential_impl<Real>(first, last); | |
441 | } | |
92f5a8d4 | 442 | } |
1e59de90 TL |
443 | |
444 | else if constexpr (std::is_integral_v<Real>) | |
445 | { | |
446 | return detail::gini_coefficient_parallel_impl<double>(exec, first, last); | |
447 | } | |
448 | ||
92f5a8d4 TL |
449 | else |
450 | { | |
1e59de90 TL |
451 | return detail::gini_coefficient_parallel_impl<Real>(exec, first, last); |
452 | } | |
453 | } | |
454 | #else | |
455 | template<class ExecutionPolicy, class RandomAccessIterator> | |
456 | inline auto gini_coefficient(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last) | |
457 | { | |
458 | using Real = typename std::iterator_traits<RandomAccessIterator>::value_type; | |
92f5a8d4 | 459 | |
1e59de90 TL |
460 | if (!std::is_sorted(exec, first, last)) |
461 | { | |
462 | std::sort(exec, first, last); | |
463 | } | |
92f5a8d4 | 464 | |
1e59de90 TL |
465 | if constexpr (std::is_integral_v<Real>) |
466 | { | |
467 | return detail::gini_coefficient_sequential_impl<double>(first, last); | |
468 | } | |
469 | else | |
470 | { | |
471 | return detail::gini_coefficient_sequential_impl<Real>(first, last); | |
472 | } | |
473 | } | |
474 | #endif | |
475 | ||
476 | template<class ExecutionPolicy, class RandomAccessContainer> | |
477 | inline auto gini_coefficient(ExecutionPolicy&& exec, RandomAccessContainer & v) | |
478 | { | |
479 | return gini_coefficient(exec, std::begin(v), std::end(v)); | |
480 | } | |
481 | ||
482 | template<class RandomAccessIterator> | |
483 | inline auto gini_coefficient(RandomAccessIterator first, RandomAccessIterator last) | |
484 | { | |
485 | return gini_coefficient(std::execution::seq, first, last); | |
92f5a8d4 TL |
486 | } |
487 | ||
488 | template<class RandomAccessContainer> | |
489 | inline auto gini_coefficient(RandomAccessContainer & v) | |
490 | { | |
1e59de90 TL |
491 | return gini_coefficient(std::execution::seq, std::begin(v), std::end(v)); |
492 | } | |
493 | ||
494 | template<class ExecutionPolicy, class RandomAccessIterator> | |
495 | inline auto sample_gini_coefficient(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last) | |
496 | { | |
497 | const auto n = std::distance(first, last); | |
498 | return n*gini_coefficient(exec, first, last)/(n-1); | |
499 | } | |
500 | ||
501 | template<class ExecutionPolicy, class RandomAccessContainer> | |
502 | inline auto sample_gini_coefficient(ExecutionPolicy&& exec, RandomAccessContainer & v) | |
503 | { | |
504 | return sample_gini_coefficient(exec, std::begin(v), std::end(v)); | |
92f5a8d4 TL |
505 | } |
506 | ||
507 | template<class RandomAccessIterator> | |
508 | inline auto sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last) | |
509 | { | |
1e59de90 | 510 | return sample_gini_coefficient(std::execution::seq, first, last); |
92f5a8d4 TL |
511 | } |
512 | ||
513 | template<class RandomAccessContainer> | |
514 | inline auto sample_gini_coefficient(RandomAccessContainer & v) | |
515 | { | |
1e59de90 | 516 | return sample_gini_coefficient(std::execution::seq, std::begin(v), std::end(v)); |
92f5a8d4 TL |
517 | } |
518 | ||
1e59de90 TL |
519 | template<class ExecutionPolicy, class RandomAccessIterator> |
520 | auto median_absolute_deviation(ExecutionPolicy&& exec, RandomAccessIterator first, RandomAccessIterator last, | |
521 | typename std::iterator_traits<RandomAccessIterator>::value_type center=std::numeric_limits<typename std::iterator_traits<RandomAccessIterator>::value_type>::quiet_NaN()) | |
92f5a8d4 TL |
522 | { |
523 | using std::abs; | |
524 | using Real = typename std::iterator_traits<RandomAccessIterator>::value_type; | |
525 | using std::isnan; | |
526 | if (isnan(center)) | |
527 | { | |
1e59de90 | 528 | center = boost::math::statistics::median(exec, first, last); |
92f5a8d4 | 529 | } |
1e59de90 TL |
530 | const auto num_elems = std::distance(first, last); |
531 | BOOST_MATH_ASSERT_MSG(num_elems > 0, "The median of a zero-length vector is undefined."); | |
92f5a8d4 TL |
532 | auto comparator = [¢er](Real a, Real b) { return abs(a-center) < abs(b-center);}; |
533 | if (num_elems & 1) | |
534 | { | |
535 | auto middle = first + (num_elems - 1)/2; | |
1e59de90 | 536 | std::nth_element(exec, first, middle, last, comparator); |
92f5a8d4 TL |
537 | return abs(*middle); |
538 | } | |
539 | else | |
540 | { | |
541 | auto middle = first + num_elems/2 - 1; | |
1e59de90 TL |
542 | std::nth_element(exec, first, middle, last, comparator); |
543 | std::nth_element(exec, middle, middle+1, last, comparator); | |
92f5a8d4 TL |
544 | return (abs(*middle) + abs(*(middle+1)))/abs(static_cast<Real>(2)); |
545 | } | |
546 | } | |
547 | ||
1e59de90 TL |
548 | template<class ExecutionPolicy, class RandomAccessContainer> |
549 | inline auto median_absolute_deviation(ExecutionPolicy&& exec, RandomAccessContainer & v, | |
550 | typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN()) | |
551 | { | |
552 | return median_absolute_deviation(exec, std::begin(v), std::end(v), center); | |
553 | } | |
554 | ||
555 | template<class RandomAccessIterator> | |
556 | inline auto median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator last, | |
557 | typename RandomAccessIterator::value_type center=std::numeric_limits<typename RandomAccessIterator::value_type>::quiet_NaN()) | |
558 | { | |
559 | return median_absolute_deviation(std::execution::seq, first, last, center); | |
560 | } | |
561 | ||
92f5a8d4 | 562 | template<class RandomAccessContainer> |
1e59de90 TL |
563 | inline auto median_absolute_deviation(RandomAccessContainer & v, |
564 | typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN()) | |
92f5a8d4 | 565 | { |
1e59de90 | 566 | return median_absolute_deviation(std::execution::seq, std::begin(v), std::end(v), center); |
92f5a8d4 TL |
567 | } |
568 | ||
1e59de90 TL |
569 | template<class ExecutionPolicy, class ForwardIterator> |
570 | auto interquartile_range(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last) | |
f67539c2 TL |
571 | { |
572 | using Real = typename std::iterator_traits<ForwardIterator>::value_type; | |
1e59de90 | 573 | static_assert(!std::is_integral_v<Real>, "Integer values have not yet been implemented."); |
f67539c2 | 574 | auto m = std::distance(first,last); |
1e59de90 | 575 | BOOST_MATH_ASSERT_MSG(m >= 3, "At least 3 samples are required to compute the interquartile range."); |
f67539c2 TL |
576 | auto k = m/4; |
577 | auto j = m - (4*k); | |
578 | // m = 4k+j. | |
579 | // If j = 0 or j = 1, then there are an even number of samples below the median, and an even number above the median. | |
580 | // Then we must average adjacent elements to get the quartiles. | |
581 | // If j = 2 or j = 3, there are an odd number of samples above and below the median, these elements may be directly extracted to get the quartiles. | |
582 | ||
583 | if (j==2 || j==3) | |
584 | { | |
585 | auto q1 = first + k; | |
586 | auto q3 = first + 3*k + j - 1; | |
1e59de90 | 587 | std::nth_element(exec, first, q1, last); |
f67539c2 | 588 | Real Q1 = *q1; |
1e59de90 | 589 | std::nth_element(exec, q1, q3, last); |
f67539c2 TL |
590 | Real Q3 = *q3; |
591 | return Q3 - Q1; | |
592 | } else { | |
593 | // j == 0 or j==1: | |
594 | auto q1 = first + k - 1; | |
595 | auto q3 = first + 3*k - 1 + j; | |
1e59de90 | 596 | std::nth_element(exec, first, q1, last); |
f67539c2 | 597 | Real a = *q1; |
1e59de90 | 598 | std::nth_element(exec, q1, q1 + 1, last); |
f67539c2 TL |
599 | Real b = *(q1 + 1); |
600 | Real Q1 = (a+b)/2; | |
1e59de90 | 601 | std::nth_element(exec, q1, q3, last); |
f67539c2 | 602 | a = *q3; |
1e59de90 | 603 | std::nth_element(exec, q3, q3 + 1, last); |
f67539c2 TL |
604 | b = *(q3 + 1); |
605 | Real Q3 = (a+b)/2; | |
606 | return Q3 - Q1; | |
607 | } | |
608 | } | |
609 | ||
1e59de90 TL |
610 | template<class ExecutionPolicy, class RandomAccessContainer> |
611 | inline auto interquartile_range(ExecutionPolicy&& exec, RandomAccessContainer & v) | |
f67539c2 | 612 | { |
1e59de90 | 613 | return interquartile_range(exec, std::begin(v), std::end(v)); |
f67539c2 TL |
614 | } |
615 | ||
1e59de90 TL |
616 | template<class RandomAccessIterator> |
617 | inline auto interquartile_range(RandomAccessIterator first, RandomAccessIterator last) | |
20effc67 | 618 | { |
1e59de90 TL |
619 | return interquartile_range(std::execution::seq, first, last); |
620 | } | |
20effc67 | 621 | |
1e59de90 TL |
622 | template<class RandomAccessContainer> |
623 | inline auto interquartile_range(RandomAccessContainer & v) | |
624 | { | |
625 | return interquartile_range(std::execution::seq, std::begin(v), std::end(v)); | |
626 | } | |
20effc67 | 627 | |
1e59de90 TL |
628 | template<class ExecutionPolicy, class ForwardIterator, class OutputIterator> |
629 | inline OutputIterator mode(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last, OutputIterator output) | |
630 | { | |
631 | if(!std::is_sorted(exec, first, last)) | |
20effc67 | 632 | { |
1e59de90 | 633 | if constexpr (std::is_same_v<typename std::iterator_traits<ForwardIterator>::iterator_category(), std::random_access_iterator_tag>) |
20effc67 | 634 | { |
1e59de90 | 635 | std::sort(exec, first, last); |
20effc67 | 636 | } |
1e59de90 | 637 | else |
20effc67 | 638 | { |
1e59de90 | 639 | BOOST_MATH_ASSERT("Data must be sorted for sequential mode calculation"); |
20effc67 | 640 | } |
1e59de90 | 641 | } |
20effc67 | 642 | |
1e59de90 TL |
643 | return detail::mode_impl(first, last, output); |
644 | } | |
645 | ||
646 | template<class ExecutionPolicy, class Container, class OutputIterator> | |
647 | inline OutputIterator mode(ExecutionPolicy&& exec, Container & v, OutputIterator output) | |
648 | { | |
649 | return mode(exec, std::begin(v), std::end(v), output); | |
650 | } | |
651 | ||
652 | template<class ForwardIterator, class OutputIterator> | |
653 | inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIterator output) | |
654 | { | |
655 | return mode(std::execution::seq, first, last, output); | |
656 | } | |
657 | ||
658 | // Requires enable_if_t to not clash with impl that returns std::list | |
659 | // Very ugly. std::is_execution_policy_v returns false for the std::execution objects and decltype of the objects (e.g. std::execution::seq) | |
660 | template<class Container, class OutputIterator, std::enable_if_t<!std::is_convertible_v<std::execution::sequenced_policy, Container> && | |
661 | !std::is_convertible_v<std::execution::parallel_unsequenced_policy, Container> && | |
662 | !std::is_convertible_v<std::execution::parallel_policy, Container> | |
663 | #if __cpp_lib_execution > 201900 | |
664 | && !std::is_convertible_v<std::execution::unsequenced_policy, Container> | |
665 | #endif | |
666 | , bool> = true> | |
667 | inline OutputIterator mode(Container & v, OutputIterator output) | |
668 | { | |
669 | return mode(std::execution::seq, std::begin(v), std::end(v), output); | |
670 | } | |
671 | ||
672 | // std::list is the return type for the proposed STL stats library | |
20effc67 | 673 | |
1e59de90 TL |
674 | template<class ExecutionPolicy, class ForwardIterator, class Real = typename std::iterator_traits<ForwardIterator>::value_type> |
675 | inline auto mode(ExecutionPolicy&& exec, ForwardIterator first, ForwardIterator last) | |
676 | { | |
677 | std::list<Real> modes; | |
678 | mode(exec, first, last, std::inserter(modes, modes.begin())); | |
679 | return modes; | |
680 | } | |
681 | ||
682 | template<class ExecutionPolicy, class Container> | |
683 | inline auto mode(ExecutionPolicy&& exec, Container & v) | |
684 | { | |
685 | return mode(exec, std::begin(v), std::end(v)); | |
686 | } | |
687 | ||
688 | template<class ForwardIterator> | |
689 | inline auto mode(ForwardIterator first, ForwardIterator last) | |
690 | { | |
691 | return mode(std::execution::seq, first, last); | |
692 | } | |
693 | ||
694 | template<class Container> | |
695 | inline auto mode(Container & v) | |
696 | { | |
697 | return mode(std::execution::seq, std::begin(v), std::end(v)); | |
698 | } | |
699 | ||
700 | } // Namespace boost::math::statistics | |
701 | ||
702 | #else // Backwards compatible bindings for C++11 | |
703 | ||
704 | namespace boost { namespace math { namespace statistics { | |
705 | ||
706 | template<bool B, class T = void> | |
707 | using enable_if_t = typename std::enable_if<B, T>::type; | |
708 | ||
709 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
710 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
711 | inline double mean(const ForwardIterator first, const ForwardIterator last) | |
712 | { | |
713 | BOOST_MATH_ASSERT_MSG(first != last, "At least one sample is required to compute the mean."); | |
714 | return detail::mean_sequential_impl<double>(first, last); | |
715 | } | |
716 | ||
717 | template<class Container, typename Real = typename Container::value_type, | |
718 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
719 | inline double mean(const Container& c) | |
720 | { | |
721 | return mean(std::begin(c), std::end(c)); | |
722 | } | |
723 | ||
724 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
725 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
726 | inline Real mean(const ForwardIterator first, const ForwardIterator last) | |
727 | { | |
728 | BOOST_MATH_ASSERT_MSG(first != last, "At least one sample is required to compute the mean."); | |
729 | return detail::mean_sequential_impl<Real>(first, last); | |
730 | } | |
731 | ||
732 | template<class Container, typename Real = typename Container::value_type, | |
733 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
734 | inline Real mean(const Container& c) | |
735 | { | |
736 | return mean(std::begin(c), std::end(c)); | |
737 | } | |
738 | ||
739 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
740 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
741 | inline double variance(const ForwardIterator first, const ForwardIterator last) | |
742 | { | |
743 | return std::get<2>(detail::variance_sequential_impl<std::tuple<double, double, double, double>>(first, last)); | |
744 | } | |
745 | ||
746 | template<class Container, typename Real = typename Container::value_type, | |
747 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
748 | inline double variance(const Container& c) | |
749 | { | |
750 | return variance(std::begin(c), std::end(c)); | |
751 | } | |
752 | ||
753 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
754 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
755 | inline Real variance(const ForwardIterator first, const ForwardIterator last) | |
756 | { | |
757 | return std::get<2>(detail::variance_sequential_impl<std::tuple<Real, Real, Real, Real>>(first, last)); | |
758 | ||
759 | } | |
760 | ||
761 | template<class Container, typename Real = typename Container::value_type, | |
762 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
763 | inline Real variance(const Container& c) | |
764 | { | |
765 | return variance(std::begin(c), std::end(c)); | |
766 | } | |
767 | ||
768 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
769 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
770 | inline double sample_variance(const ForwardIterator first, const ForwardIterator last) | |
771 | { | |
772 | const auto n = std::distance(first, last); | |
773 | BOOST_MATH_ASSERT_MSG(n > 1, "At least two samples are required to compute the sample variance."); | |
774 | return n*variance(first, last)/(n-1); | |
775 | } | |
776 | ||
777 | template<class Container, typename Real = typename Container::value_type, | |
778 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
779 | inline double sample_variance(const Container& c) | |
780 | { | |
781 | return sample_variance(std::begin(c), std::end(c)); | |
782 | } | |
783 | ||
784 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
785 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
786 | inline Real sample_variance(const ForwardIterator first, const ForwardIterator last) | |
787 | { | |
788 | const auto n = std::distance(first, last); | |
789 | BOOST_MATH_ASSERT_MSG(n > 1, "At least two samples are required to compute the sample variance."); | |
790 | return n*variance(first, last)/(n-1); | |
791 | } | |
792 | ||
793 | template<class Container, typename Real = typename Container::value_type, | |
794 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
795 | inline Real sample_variance(const Container& c) | |
796 | { | |
797 | return sample_variance(std::begin(c), std::end(c)); | |
798 | } | |
799 | ||
800 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
801 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
802 | inline std::pair<double, double> mean_and_sample_variance(const ForwardIterator first, const ForwardIterator last) | |
803 | { | |
804 | const auto results = detail::variance_sequential_impl<std::tuple<double, double, double, double>>(first, last); | |
805 | return std::make_pair(std::get<0>(results), std::get<3>(results)*std::get<2>(results)/(std::get<3>(results)-1.0)); | |
806 | } | |
807 | ||
808 | template<class Container, typename Real = typename Container::value_type, | |
809 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
810 | inline std::pair<double, double> mean_and_sample_variance(const Container& c) | |
811 | { | |
812 | return mean_and_sample_variance(std::begin(c), std::end(c)); | |
813 | } | |
814 | ||
815 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
816 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
817 | inline std::pair<Real, Real> mean_and_sample_variance(const ForwardIterator first, const ForwardIterator last) | |
818 | { | |
819 | const auto results = detail::variance_sequential_impl<std::tuple<Real, Real, Real, Real>>(first, last); | |
820 | return std::make_pair(std::get<0>(results), std::get<3>(results)*std::get<2>(results)/(std::get<3>(results)-Real(1))); | |
821 | } | |
822 | ||
823 | template<class Container, typename Real = typename Container::value_type, | |
824 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
825 | inline std::pair<Real, Real> mean_and_sample_variance(const Container& c) | |
826 | { | |
827 | return mean_and_sample_variance(std::begin(c), std::end(c)); | |
828 | } | |
829 | ||
830 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
831 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
832 | inline std::tuple<double, double, double, double> first_four_moments(const ForwardIterator first, const ForwardIterator last) | |
833 | { | |
834 | const auto results = detail::first_four_moments_sequential_impl<std::tuple<double, double, double, double, double>>(first, last); | |
835 | return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results), | |
836 | std::get<3>(results) / std::get<4>(results)); | |
837 | } | |
838 | ||
839 | template<class Container, typename Real = typename Container::value_type, | |
840 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
841 | inline std::tuple<double, double, double, double> first_four_moments(const Container& c) | |
842 | { | |
843 | return first_four_moments(std::begin(c), std::end(c)); | |
844 | } | |
845 | ||
846 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
847 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
848 | inline std::tuple<Real, Real, Real, Real> first_four_moments(const ForwardIterator first, const ForwardIterator last) | |
849 | { | |
850 | const auto results = detail::first_four_moments_sequential_impl<std::tuple<Real, Real, Real, Real, Real>>(first, last); | |
851 | return std::make_tuple(std::get<0>(results), std::get<1>(results) / std::get<4>(results), std::get<2>(results) / std::get<4>(results), | |
852 | std::get<3>(results) / std::get<4>(results)); | |
853 | } | |
854 | ||
855 | template<class Container, typename Real = typename Container::value_type, | |
856 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
857 | inline std::tuple<Real, Real, Real, Real> first_four_moments(const Container& c) | |
858 | { | |
859 | return first_four_moments(std::begin(c), std::end(c)); | |
860 | } | |
861 | ||
862 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
863 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
864 | inline double skewness(const ForwardIterator first, const ForwardIterator last) | |
865 | { | |
866 | return detail::skewness_sequential_impl<double>(first, last); | |
867 | } | |
868 | ||
869 | template<class Container, typename Real = typename Container::value_type, | |
870 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
871 | inline double skewness(const Container& c) | |
872 | { | |
873 | return skewness(std::begin(c), std::end(c)); | |
874 | } | |
875 | ||
876 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
877 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
878 | inline Real skewness(const ForwardIterator first, const ForwardIterator last) | |
879 | { | |
880 | return detail::skewness_sequential_impl<Real>(first, last); | |
881 | } | |
882 | ||
883 | template<class Container, typename Real = typename Container::value_type, | |
884 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
885 | inline Real skewness(const Container& c) | |
886 | { | |
887 | return skewness(std::begin(c), std::end(c)); | |
888 | } | |
889 | ||
890 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
891 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
892 | inline double kurtosis(const ForwardIterator first, const ForwardIterator last) | |
893 | { | |
894 | std::tuple<double, double, double, double> M = first_four_moments(first, last); | |
895 | ||
896 | if(std::get<1>(M) == 0) | |
897 | { | |
898 | return std::get<1>(M); | |
899 | } | |
900 | else | |
901 | { | |
902 | return std::get<3>(M)/(std::get<1>(M)*std::get<1>(M)); | |
20effc67 | 903 | } |
1e59de90 | 904 | } |
20effc67 | 905 | |
1e59de90 TL |
906 | template<class Container, typename Real = typename Container::value_type, |
907 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
908 | inline double kurtosis(const Container& c) | |
909 | { | |
910 | return kurtosis(std::begin(c), std::end(c)); | |
20effc67 TL |
911 | } |
912 | ||
1e59de90 TL |
913 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, |
914 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
915 | inline Real kurtosis(const ForwardIterator first, const ForwardIterator last) | |
916 | { | |
917 | std::tuple<Real, Real, Real, Real> M = first_four_moments(first, last); | |
918 | ||
919 | if(std::get<1>(M) == 0) | |
920 | { | |
921 | return std::get<1>(M); | |
922 | } | |
923 | else | |
924 | { | |
925 | return std::get<3>(M)/(std::get<1>(M)*std::get<1>(M)); | |
926 | } | |
927 | } | |
928 | ||
929 | template<class Container, typename Real = typename Container::value_type, | |
930 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
931 | inline Real kurtosis(const Container& c) | |
932 | { | |
933 | return kurtosis(std::begin(c), std::end(c)); | |
934 | } | |
935 | ||
936 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
937 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
938 | inline double excess_kurtosis(const ForwardIterator first, const ForwardIterator last) | |
20effc67 | 939 | { |
1e59de90 TL |
940 | return kurtosis(first, last) - 3; |
941 | } | |
942 | ||
943 | template<class Container, typename Real = typename Container::value_type, | |
944 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
945 | inline double excess_kurtosis(const Container& c) | |
946 | { | |
947 | return excess_kurtosis(std::begin(c), std::end(c)); | |
948 | } | |
949 | ||
950 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type, | |
951 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
952 | inline Real excess_kurtosis(const ForwardIterator first, const ForwardIterator last) | |
953 | { | |
954 | return kurtosis(first, last) - 3; | |
955 | } | |
956 | ||
957 | template<class Container, typename Real = typename Container::value_type, | |
958 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
959 | inline Real excess_kurtosis(const Container& c) | |
960 | { | |
961 | return excess_kurtosis(std::begin(c), std::end(c)); | |
20effc67 TL |
962 | } |
963 | ||
1e59de90 TL |
964 | template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type> |
965 | Real median(RandomAccessIterator first, RandomAccessIterator last) | |
20effc67 | 966 | { |
1e59de90 TL |
967 | const auto num_elems = std::distance(first, last); |
968 | BOOST_MATH_ASSERT_MSG(num_elems > 0, "The median of a zero length vector is undefined."); | |
969 | if (num_elems & 1) | |
970 | { | |
971 | auto middle = first + (num_elems - 1)/2; | |
972 | std::nth_element(first, middle, last); | |
973 | return *middle; | |
974 | } | |
975 | else | |
976 | { | |
977 | auto middle = first + num_elems/2 - 1; | |
978 | std::nth_element(first, middle, last); | |
979 | std::nth_element(middle, middle+1, last); | |
980 | return (*middle + *(middle+1))/2; | |
981 | } | |
20effc67 TL |
982 | } |
983 | ||
1e59de90 TL |
984 | template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type> |
985 | inline Real median(RandomAccessContainer& c) | |
20effc67 | 986 | { |
1e59de90 | 987 | return median(std::begin(c), std::end(c)); |
20effc67 | 988 | } |
f67539c2 | 989 | |
1e59de90 TL |
990 | template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type, |
991 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
992 | inline double gini_coefficient(RandomAccessIterator first, RandomAccessIterator last) | |
993 | { | |
994 | if(!std::is_sorted(first, last)) | |
995 | { | |
996 | std::sort(first, last); | |
997 | } | |
998 | ||
999 | return detail::gini_coefficient_sequential_impl<double>(first, last); | |
1000 | } | |
1001 | ||
1002 | template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type, | |
1003 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
1004 | inline double gini_coefficient(RandomAccessContainer& c) | |
1005 | { | |
1006 | return gini_coefficient(std::begin(c), std::end(c)); | |
1007 | } | |
1008 | ||
1009 | template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type, | |
1010 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
1011 | inline Real gini_coefficient(RandomAccessIterator first, RandomAccessIterator last) | |
1012 | { | |
1013 | if(!std::is_sorted(first, last)) | |
1014 | { | |
1015 | std::sort(first, last); | |
1016 | } | |
1017 | ||
1018 | return detail::gini_coefficient_sequential_impl<Real>(first, last); | |
1019 | } | |
1020 | ||
1021 | template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type, | |
1022 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
1023 | inline Real gini_coefficient(RandomAccessContainer& c) | |
1024 | { | |
1025 | return gini_coefficient(std::begin(c), std::end(c)); | |
1026 | } | |
1027 | ||
1028 | template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type, | |
1029 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
1030 | inline double sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last) | |
1031 | { | |
1032 | const auto n = std::distance(first, last); | |
1033 | return n*gini_coefficient(first, last)/(n-1); | |
1034 | } | |
1035 | ||
1036 | template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type, | |
1037 | enable_if_t<std::is_integral<Real>::value, bool> = true> | |
1038 | inline double sample_gini_coefficient(RandomAccessContainer& c) | |
1039 | { | |
1040 | return sample_gini_coefficient(std::begin(c), std::end(c)); | |
1041 | } | |
1042 | ||
1043 | template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type, | |
1044 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
1045 | inline Real sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last) | |
1046 | { | |
1047 | const auto n = std::distance(first, last); | |
1048 | return n*gini_coefficient(first, last)/(n-1); | |
1049 | } | |
1050 | ||
1051 | template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type, | |
1052 | enable_if_t<!std::is_integral<Real>::value, bool> = true> | |
1053 | inline Real sample_gini_coefficient(RandomAccessContainer& c) | |
1054 | { | |
1055 | return sample_gini_coefficient(std::begin(c), std::end(c)); | |
1056 | } | |
1057 | ||
1058 | template<class RandomAccessIterator, typename Real = typename std::iterator_traits<RandomAccessIterator>::value_type> | |
1059 | Real median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator last, | |
1060 | typename std::iterator_traits<RandomAccessIterator>::value_type center=std::numeric_limits<typename std::iterator_traits<RandomAccessIterator>::value_type>::quiet_NaN()) | |
1061 | { | |
1062 | using std::abs; | |
1063 | using std::isnan; | |
1064 | if (isnan(center)) | |
1065 | { | |
1066 | center = boost::math::statistics::median(first, last); | |
1067 | } | |
1068 | const auto num_elems = std::distance(first, last); | |
1069 | BOOST_MATH_ASSERT_MSG(num_elems > 0, "The median of a zero-length vector is undefined."); | |
1070 | auto comparator = [¢er](Real a, Real b) { return abs(a-center) < abs(b-center);}; | |
1071 | if (num_elems & 1) | |
1072 | { | |
1073 | auto middle = first + (num_elems - 1)/2; | |
1074 | std::nth_element(first, middle, last, comparator); | |
1075 | return abs(*middle); | |
1076 | } | |
1077 | else | |
1078 | { | |
1079 | auto middle = first + num_elems/2 - 1; | |
1080 | std::nth_element(first, middle, last, comparator); | |
1081 | std::nth_element(middle, middle+1, last, comparator); | |
1082 | return (abs(*middle) + abs(*(middle+1)))/abs(static_cast<Real>(2)); | |
1083 | } | |
1084 | } | |
1085 | ||
1086 | template<class RandomAccessContainer, typename Real = typename RandomAccessContainer::value_type> | |
1087 | inline Real median_absolute_deviation(RandomAccessContainer& c, | |
1088 | typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN()) | |
1089 | { | |
1090 | return median_absolute_deviation(std::begin(c), std::end(c), center); | |
1091 | } | |
1092 | ||
1093 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type> | |
1094 | Real interquartile_range(ForwardIterator first, ForwardIterator last) | |
1095 | { | |
1096 | static_assert(!std::is_integral<Real>::value, "Integer values have not yet been implemented."); | |
1097 | auto m = std::distance(first,last); | |
1098 | BOOST_MATH_ASSERT_MSG(m >= 3, "At least 3 samples are required to compute the interquartile range."); | |
1099 | auto k = m/4; | |
1100 | auto j = m - (4*k); | |
1101 | // m = 4k+j. | |
1102 | // If j = 0 or j = 1, then there are an even number of samples below the median, and an even number above the median. | |
1103 | // Then we must average adjacent elements to get the quartiles. | |
1104 | // If j = 2 or j = 3, there are an odd number of samples above and below the median, these elements may be directly extracted to get the quartiles. | |
1105 | ||
1106 | if (j==2 || j==3) | |
1107 | { | |
1108 | auto q1 = first + k; | |
1109 | auto q3 = first + 3*k + j - 1; | |
1110 | std::nth_element(first, q1, last); | |
1111 | Real Q1 = *q1; | |
1112 | std::nth_element(q1, q3, last); | |
1113 | Real Q3 = *q3; | |
1114 | return Q3 - Q1; | |
1115 | } | |
1116 | else | |
1117 | { | |
1118 | // j == 0 or j==1: | |
1119 | auto q1 = first + k - 1; | |
1120 | auto q3 = first + 3*k - 1 + j; | |
1121 | std::nth_element(first, q1, last); | |
1122 | Real a = *q1; | |
1123 | std::nth_element(q1, q1 + 1, last); | |
1124 | Real b = *(q1 + 1); | |
1125 | Real Q1 = (a+b)/2; | |
1126 | std::nth_element(q1, q3, last); | |
1127 | a = *q3; | |
1128 | std::nth_element(q3, q3 + 1, last); | |
1129 | b = *(q3 + 1); | |
1130 | Real Q3 = (a+b)/2; | |
1131 | return Q3 - Q1; | |
1132 | } | |
1133 | } | |
1134 | ||
1135 | template<class Container, typename Real = typename Container::value_type> | |
1136 | Real interquartile_range(Container& c) | |
1137 | { | |
1138 | return interquartile_range(std::begin(c), std::end(c)); | |
1139 | } | |
1140 | ||
1141 | template<class ForwardIterator, class OutputIterator, | |
1142 | enable_if_t<std::is_same<typename std::iterator_traits<ForwardIterator>::iterator_category(), std::random_access_iterator_tag>::value, bool> = true> | |
1143 | inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIterator output) | |
1144 | { | |
1145 | if(!std::is_sorted(first, last)) | |
1146 | { | |
1147 | std::sort(first, last); | |
1148 | } | |
1149 | ||
1150 | return detail::mode_impl(first, last, output); | |
1151 | } | |
1152 | ||
1153 | template<class ForwardIterator, class OutputIterator, | |
1154 | enable_if_t<!std::is_same<typename std::iterator_traits<ForwardIterator>::iterator_category(), std::random_access_iterator_tag>::value, bool> = true> | |
1155 | inline OutputIterator mode(ForwardIterator first, ForwardIterator last, OutputIterator output) | |
1156 | { | |
1157 | if(!std::is_sorted(first, last)) | |
1158 | { | |
1159 | BOOST_MATH_ASSERT("Data must be sorted for mode calculation"); | |
1160 | } | |
1161 | ||
1162 | return detail::mode_impl(first, last, output); | |
1163 | } | |
1164 | ||
1165 | template<class Container, class OutputIterator> | |
1166 | inline OutputIterator mode(Container& c, OutputIterator output) | |
1167 | { | |
1168 | return mode(std::begin(c), std::end(c), output); | |
1169 | } | |
1170 | ||
1171 | template<class ForwardIterator, typename Real = typename std::iterator_traits<ForwardIterator>::value_type> | |
1172 | inline std::list<Real> mode(ForwardIterator first, ForwardIterator last) | |
1173 | { | |
1174 | std::list<Real> modes; | |
1175 | mode(first, last, std::inserter(modes, modes.begin())); | |
1176 | return modes; | |
1177 | } | |
1178 | ||
1179 | template<class Container, typename Real = typename Container::value_type> | |
1180 | inline std::list<Real> mode(Container& c) | |
1181 | { | |
1182 | return mode(std::begin(c), std::end(c)); | |
92f5a8d4 | 1183 | } |
1e59de90 | 1184 | }}} |
92f5a8d4 | 1185 | #endif |
1e59de90 | 1186 | #endif // BOOST_MATH_STATISTICS_UNIVARIATE_STATISTICS_HPP |