boost/math/statistics/univariate_statistics.hpp

   1 //  (C) Copyright Nick Thompson 2018.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_STATISTICS_UNIVARIATE_STATISTICS_HPP
   7 #define BOOST_MATH_STATISTICS_UNIVARIATE_STATISTICS_HPP
   8
   9 #include <algorithm>
  10 #include <iterator>
  11 #include <boost/assert.hpp>
  12
  13 namespace boost::math::statistics {
  14
  15 template<class ForwardIterator>
  16 auto mean(ForwardIterator first, ForwardIterator last)
  17 {
  18     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
  19     BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute the mean.");
  20     if constexpr (std::is_integral<Real>::value)
  21     {
  22         double mu = 0;
  23         double i = 1;
  24         for(auto it = first; it != last; ++it) {
  25             mu = mu + (*it - mu)/i;
  26             i += 1;
  27         }
  28         return mu;
  29     }
  30     else if constexpr (std::is_same_v<typename std::iterator_traits<ForwardIterator>::iterator_category, std::random_access_iterator_tag>)
  31     {
  32         size_t elements = std::distance(first, last);
  33         Real mu0 = 0;
  34         Real mu1 = 0;
  35         Real mu2 = 0;
  36         Real mu3 = 0;
  37         Real i = 1;
  38         auto end = last - (elements % 4);
  39         for(auto it = first; it != end;  it += 4) {
  40             Real inv = Real(1)/i;
  41             Real tmp0 = (*it - mu0);
  42             Real tmp1 = (*(it+1) - mu1);
  43             Real tmp2 = (*(it+2) - mu2);
  44             Real tmp3 = (*(it+3) - mu3);
  45             // please generate a vectorized fma here
  46             mu0 += tmp0*inv;
  47             mu1 += tmp1*inv;
  48             mu2 += tmp2*inv;
  49             mu3 += tmp3*inv;
  50             i += 1;
  51         }
  52         Real num1 = Real(elements  - (elements %4))/Real(4);
  53         Real num2 = num1 + Real(elements % 4);
  54
  55         for (auto it = end; it != last; ++it)
  56         {
  57             mu3 += (*it-mu3)/i;
  58             i += 1;
  59         }
  60
  61         return (num1*(mu0+mu1+mu2) + num2*mu3)/Real(elements);
  62     }
  63     else
  64     {
  65         auto it = first;
  66         Real mu = *it;
  67         Real i = 2;
  68         while(++it != last)
  69         {
  70             mu += (*it - mu)/i;
  71             i += 1;
  72         }
  73         return mu;
  74     }
  75 }
  76
  77 template<class Container>
  78 inline auto mean(Container const & v)
  79 {
  80     return mean(v.cbegin(), v.cend());
  81 }
  82
  83 template<class ForwardIterator>
  84 auto variance(ForwardIterator first, ForwardIterator last)
  85 {
  86     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
  87     BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute mean and variance.");
  88     // Higham, Accuracy and Stability, equation 1.6a and 1.6b:
  89     if constexpr (std::is_integral<Real>::value)
  90     {
  91         double M = *first;
  92         double Q = 0;
  93         double k = 2;
  94         for (auto it = std::next(first); it != last; ++it)
  95         {
  96             double tmp = *it - M;
  97             Q = Q + ((k-1)*tmp*tmp)/k;
  98             M = M + tmp/k;
  99             k += 1;
 100         }
 101         return Q/(k-1);
 102     }
 103     else
 104     {
 105         Real M = *first;
 106         Real Q = 0;
 107         Real k = 2;
 108         for (auto it = std::next(first); it != last; ++it)
 109         {
 110             Real tmp = (*it - M)/k;
 111             Q += k*(k-1)*tmp*tmp;
 112             M += tmp;
 113             k += 1;
 114         }
 115         return Q/(k-1);
 116     }
 117 }
 118
 119 template<class Container>
 120 inline auto variance(Container const & v)
 121 {
 122     return variance(v.cbegin(), v.cend());
 123 }
 124
 125 template<class ForwardIterator>
 126 auto sample_variance(ForwardIterator first, ForwardIterator last)
 127 {
 128     size_t n = std::distance(first, last);
 129     BOOST_ASSERT_MSG(n > 1, "At least two samples are required to compute the sample variance.");
 130     return n*variance(first, last)/(n-1);
 131 }
 132
 133 template<class Container>
 134 inline auto sample_variance(Container const & v)
 135 {
 136     return sample_variance(v.cbegin(), v.cend());
 137 }
 138
 139 template<class ForwardIterator>
 140 auto mean_and_sample_variance(ForwardIterator first, ForwardIterator last)
 141 {
 142     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
 143     BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute mean and variance.");
 144     // Higham, Accuracy and Stability, equation 1.6a and 1.6b:
 145     if constexpr (std::is_integral<Real>::value)
 146     {
 147         double M = *first;
 148         double Q = 0;
 149         double k = 2;
 150         for (auto it = std::next(first); it != last; ++it)
 151         {
 152             double tmp = *it - M;
 153             Q = Q + ((k-1)*tmp*tmp)/k;
 154             M = M + tmp/k;
 155             k += 1;
 156         }
 157         return std::pair<double, double>{M, Q/(k-2)};
 158     }
 159     else
 160     {
 161         Real M = *first;
 162         Real Q = 0;
 163         Real k = 2;
 164         for (auto it = std::next(first); it != last; ++it)
 165         {
 166             Real tmp = (*it - M)/k;
 167             Q += k*(k-1)*tmp*tmp;
 168             M += tmp;
 169             k += 1;
 170         }
 171         return std::pair<Real, Real>{M, Q/(k-2)};
 172     }
 173 }
 174
 175 template<class Container>
 176 auto mean_and_sample_variance(Container const & v)
 177 {
 178     return mean_and_sample_variance(v.begin(), v.end());
 179 }
 180
 181 // Follows equation 1.5 of:
 182 // https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
 183 template<class ForwardIterator>
 184 auto skewness(ForwardIterator first, ForwardIterator last)
 185 {
 186     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
 187     BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute skewness.");
 188     if constexpr (std::is_integral<Real>::value)
 189     {
 190         double M1 = *first;
 191         double M2 = 0;
 192         double M3 = 0;
 193         double n = 2;
 194         for (auto it = std::next(first); it != last; ++it)
 195         {
 196             double delta21 = *it - M1;
 197             double tmp = delta21/n;
 198             M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
 199             M2 = M2 + tmp*(n-1)*delta21;
 200             M1 = M1 + tmp;
 201             n += 1;
 202         }
 203
 204         double var = M2/(n-1);
 205         if (var == 0)
 206         {
 207             // The limit is technically undefined, but the interpretation here is clear:
 208             // A constant dataset has no skewness.
 209             return double(0);
 210         }
 211         double skew = M3/(M2*sqrt(var));
 212         return skew;
 213     }
 214     else
 215     {
 216         Real M1 = *first;
 217         Real M2 = 0;
 218         Real M3 = 0;
 219         Real n = 2;
 220         for (auto it = std::next(first); it != last; ++it)
 221         {
 222             Real delta21 = *it - M1;
 223             Real tmp = delta21/n;
 224             M3 += tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
 225             M2 += tmp*(n-1)*delta21;
 226             M1 += tmp;
 227             n += 1;
 228         }
 229
 230         Real var = M2/(n-1);
 231         if (var == 0)
 232         {
 233             // The limit is technically undefined, but the interpretation here is clear:
 234             // A constant dataset has no skewness.
 235             return Real(0);
 236         }
 237         Real skew = M3/(M2*sqrt(var));
 238         return skew;
 239     }
 240 }
 241
 242 template<class Container>
 243 inline auto skewness(Container const & v)
 244 {
 245     return skewness(v.cbegin(), v.cend());
 246 }
 247
 248 // Follows equation 1.5/1.6 of:
 249 // https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
 250 template<class ForwardIterator>
 251 auto first_four_moments(ForwardIterator first, ForwardIterator last)
 252 {
 253     using Real = typename std::iterator_traits<ForwardIterator>::value_type;
 254     BOOST_ASSERT_MSG(first != last, "At least one sample is required to compute the first four moments.");
 255     if constexpr (std::is_integral<Real>::value)
 256     {
 257         double M1 = *first;
 258         double M2 = 0;
 259         double M3 = 0;
 260         double M4 = 0;
 261         double n = 2;
 262         for (auto it = std::next(first); it != last; ++it)
 263         {
 264             double delta21 = *it - M1;
 265             double tmp = delta21/n;
 266             M4 = M4 + tmp*(tmp*tmp*delta21*((n-1)*(n*n-3*n+3)) + 6*tmp*M2 - 4*M3);
 267             M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
 268             M2 = M2 + tmp*(n-1)*delta21;
 269             M1 = M1 + tmp;
 270             n += 1;
 271         }
 272
 273         return std::make_tuple(M1, M2/(n-1), M3/(n-1), M4/(n-1));
 274     }
 275     else
 276     {
 277         Real M1 = *first;
 278         Real M2 = 0;
 279         Real M3 = 0;
 280         Real M4 = 0;
 281         Real n = 2;
 282         for (auto it = std::next(first); it != last; ++it)
 283         {
 284             Real delta21 = *it - M1;
 285             Real tmp = delta21/n;
 286             M4 = M4 + tmp*(tmp*tmp*delta21*((n-1)*(n*n-3*n+3)) + 6*tmp*M2 - 4*M3);
 287             M3 = M3 + tmp*((n-1)*(n-2)*delta21*tmp - 3*M2);
 288             M2 = M2 + tmp*(n-1)*delta21;
 289             M1 = M1 + tmp;
 290             n += 1;
 291         }
 292
 293         return std::make_tuple(M1, M2/(n-1), M3/(n-1), M4/(n-1));
 294     }
 295 }
 296
 297 template<class Container>
 298 inline auto first_four_moments(Container const & v)
 299 {
 300     return first_four_moments(v.cbegin(), v.cend());
 301 }
 302
 303
 304 // Follows equation 1.6 of:
 305 // https://prod.sandia.gov/techlib-noauth/access-control.cgi/2008/086212.pdf
 306 template<class ForwardIterator>
 307 auto kurtosis(ForwardIterator first, ForwardIterator last)
 308 {
 309     auto [M1, M2, M3, M4] = first_four_moments(first, last);
 310     if (M2 == 0)
 311     {
 312         return M2;
 313     }
 314     return M4/(M2*M2);
 315 }
 316
 317 template<class Container>
 318 inline auto kurtosis(Container const & v)
 319 {
 320     return kurtosis(v.cbegin(), v.cend());
 321 }
 322
 323 template<class ForwardIterator>
 324 auto excess_kurtosis(ForwardIterator first, ForwardIterator last)
 325 {
 326     return kurtosis(first, last) - 3;
 327 }
 328
 329 template<class Container>
 330 inline auto excess_kurtosis(Container const & v)
 331 {
 332     return excess_kurtosis(v.cbegin(), v.cend());
 333 }
 334
 335
 336 template<class RandomAccessIterator>
 337 auto median(RandomAccessIterator first, RandomAccessIterator last)
 338 {
 339     size_t num_elems = std::distance(first, last);
 340     BOOST_ASSERT_MSG(num_elems > 0, "The median of a zero length vector is undefined.");
 341     if (num_elems & 1)
 342     {
 343         auto middle = first + (num_elems - 1)/2;
 344         std::nth_element(first, middle, last);
 345         return *middle;
 346     }
 347     else
 348     {
 349         auto middle = first + num_elems/2 - 1;
 350         std::nth_element(first, middle, last);
 351         std::nth_element(middle, middle+1, last);
 352         return (*middle + *(middle+1))/2;
 353     }
 354 }
 355
 356
 357 template<class RandomAccessContainer>
 358 inline auto median(RandomAccessContainer & v)
 359 {
 360     return median(v.begin(), v.end());
 361 }
 362
 363 template<class RandomAccessIterator>
 364 auto gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
 365 {
 366     using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;
 367     BOOST_ASSERT_MSG(first != last && std::next(first) != last, "Computation of the Gini coefficient requires at least two samples.");
 368
 369     std::sort(first, last);
 370     if constexpr (std::is_integral<Real>::value)
 371     {
 372         double i = 1;
 373         double num = 0;
 374         double denom = 0;
 375         for (auto it = first; it != last; ++it)
 376         {
 377             num += *it*i;
 378             denom += *it;
 379             ++i;
 380         }
 381
 382         // If the l1 norm is zero, all elements are zero, so every element is the same.
 383         if (denom == 0)
 384         {
 385             return double(0);
 386         }
 387
 388         return ((2*num)/denom - i)/(i-1);
 389     }
 390     else
 391     {
 392         Real i = 1;
 393         Real num = 0;
 394         Real denom = 0;
 395         for (auto it = first; it != last; ++it)
 396         {
 397             num += *it*i;
 398             denom += *it;
 399             ++i;
 400         }
 401
 402         // If the l1 norm is zero, all elements are zero, so every element is the same.
 403         if (denom == 0)
 404         {
 405             return Real(0);
 406         }
 407
 408         return ((2*num)/denom - i)/(i-1);
 409     }
 410 }
 411
 412 template<class RandomAccessContainer>
 413 inline auto gini_coefficient(RandomAccessContainer & v)
 414 {
 415     return gini_coefficient(v.begin(), v.end());
 416 }
 417
 418 template<class RandomAccessIterator>
 419 inline auto sample_gini_coefficient(RandomAccessIterator first, RandomAccessIterator last)
 420 {
 421     size_t n = std::distance(first, last);
 422     return n*gini_coefficient(first, last)/(n-1);
 423 }
 424
 425 template<class RandomAccessContainer>
 426 inline auto sample_gini_coefficient(RandomAccessContainer & v)
 427 {
 428     return sample_gini_coefficient(v.begin(), v.end());
 429 }
 430
 431 template<class RandomAccessIterator>
 432 auto median_absolute_deviation(RandomAccessIterator first, RandomAccessIterator last, typename std::iterator_traits<RandomAccessIterator>::value_type center=std::numeric_limits<typename std::iterator_traits<RandomAccessIterator>::value_type>::quiet_NaN())
 433 {
 434     using std::abs;
 435     using Real = typename std::iterator_traits<RandomAccessIterator>::value_type;
 436     using std::isnan;
 437     if (isnan(center))
 438     {
 439         center = boost::math::statistics::median(first, last);
 440     }
 441     size_t num_elems = std::distance(first, last);
 442     BOOST_ASSERT_MSG(num_elems > 0, "The median of a zero-length vector is undefined.");
 443     auto comparator = [&center](Real a, Real b) { return abs(a-center) < abs(b-center);};
 444     if (num_elems & 1)
 445     {
 446         auto middle = first + (num_elems - 1)/2;
 447         std::nth_element(first, middle, last, comparator);
 448         return abs(*middle);
 449     }
 450     else
 451     {
 452         auto middle = first + num_elems/2 - 1;
 453         std::nth_element(first, middle, last, comparator);
 454         std::nth_element(middle, middle+1, last, comparator);
 455         return (abs(*middle) + abs(*(middle+1)))/abs(static_cast<Real>(2));
 456     }
 457 }
 458
 459 template<class RandomAccessContainer>
 460 inline auto median_absolute_deviation(RandomAccessContainer & v, typename RandomAccessContainer::value_type center=std::numeric_limits<typename RandomAccessContainer::value_type>::quiet_NaN())
 461 {
 462     return median_absolute_deviation(v.begin(), v.end(), center);
 463 }
 464
 465 }
 466 #endif