StatisticFunctions.h

Go to the documentation of this file.

// Copyright (c) 2002-present, OpenMS Inc. -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg $
// $Authors: Clemens Groepl, Johannes Junker, Mathias Walzer, Chris Bielow $
// --------------------------------------------------------------------------
#pragma once
 
#include <vector>
#include <OpenMS/CONCEPT/Exception.h>
#include <OpenMS/CONCEPT/Macros.h>
#include <OpenMS/CONCEPT/Types.h>
#include <OpenMS/DATASTRUCTURES/StringUtils.h>
 
#include <algorithm>
#include <cmath>
#include <iterator>
#include <numeric>
 
namespace OpenMS
{
 
  namespace Math
  {
    struct AdaptiveQuantileResult
    {
      double blended{0.0};
      double half_raw{0.0};
      double half_rob{0.0};
      double upper_fence{std::numeric_limits<double>::infinity()};
      double tail_fraction{0.0};
      double weight{0.0};
    };
 
    template <typename IteratorType>
    static void checkIteratorsNotNULL(IteratorType begin, IteratorType end)
    {
      if (begin == end)
      {
        throw Exception::InvalidRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION);
      }
    }
 
    template <typename IteratorType>
    static void checkIteratorsEqual(IteratorType begin, IteratorType end)
    {
      if (begin != end)
      {
        throw Exception::InvalidRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION);
      }
    }
 
    template <typename IteratorType1, typename IteratorType2>
    static void checkIteratorsAreValid(
      IteratorType1 begin_b, IteratorType1 end_b,
      IteratorType2 begin_a, IteratorType2 end_a)
    {
      if ((begin_b == end_b) ^ (begin_a == end_a))
      {
        throw Exception::InvalidRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION);
      }
    }
 
    template <typename IteratorType>
    static double sum(IteratorType begin, IteratorType end)
    {
      return std::accumulate(begin, end, 0.0);
    }
 
    template <typename IteratorType>
    static double mean(IteratorType begin, IteratorType end)
    {
      checkIteratorsNotNULL(begin, end);
      return sum(begin, end) / std::distance(begin, end);
    }
 
    template <typename IteratorType>
    static double median(IteratorType begin, IteratorType end, 
                         bool sorted = false)
    {
      checkIteratorsNotNULL(begin, end);
      if (!sorted)
      {
        std::sort(begin, end);
      }
      
      Size size = std::distance(begin, end);
      if (size % 2 == 0) // even size => average two middle values
      {
        IteratorType it1 = begin;
        std::advance(it1, size / 2 - 1);
        IteratorType it2 = it1;
        std::advance(it2, 1);
        return (*it1 + *it2) / 2.0;
      }
      else
      {
        IteratorType it = begin;
        std::advance(it, (size - 1) / 2);
        return *it;
      }
    }
 
  
    template <typename IteratorType>
    double MAD(IteratorType begin, IteratorType end, double median_of_numbers)
    {
      std::vector<double> diffs;
      diffs.reserve(std::distance(begin, end));
      for (IteratorType it = begin; it != end; ++it)
      {
        diffs.push_back(fabs(*it - median_of_numbers));
      }
      return median(diffs.begin(), diffs.end(), false);
    }
    
    template <typename IteratorType>
    double MeanAbsoluteDeviation(IteratorType begin, IteratorType end, double mean_of_numbers)
    {
      double mean_value {0};
      for (IteratorType it = begin; it != end; ++it)
      {
        mean_value += fabs(*it - mean_of_numbers);
      }
      return mean_value / std::distance(begin, end);
    }
 
    template <typename IteratorType>
    static double quantile1st(IteratorType begin, IteratorType end,
                              bool sorted = false)
    {
      checkIteratorsNotNULL(begin, end);
 
      if (!sorted)
      {
        std::sort(begin, end);
      }
 
      Size size = std::distance(begin, end);
      if (size < 3) // lower half is empty for size 1 or 2: the first quantile is the minimum
      {
        return static_cast<double>(*begin);
      }
      if (size % 2 == 0)
      {
        return median(begin, begin + (size/2)-1, true); //-1 to exclude median values
      }
      return median(begin, begin + (size/2), true);
    }
 
    template <typename IteratorType>
    static double quantile3rd(
      IteratorType begin, IteratorType end, bool sorted = false)
    {
      checkIteratorsNotNULL(begin, end);
      if (!sorted)
      {
        std::sort(begin, end);
      }
 
      Size size = std::distance(begin, end);
      if (size < 3) // upper half is empty for size 1 or 2: the third quantile is the maximum
      {
        return static_cast<double>(*(begin + (size - 1)));
      }
      return median(begin + (size/2)+1, end, true); //+1 to exclude median values
    }
 
    template <typename IteratorType>
    static double quantile(IteratorType begin, IteratorType end, double q)
    {
      OPENMS_PRECONDITION(std::is_sorted(begin, end),
                          "Math::quantile expects a sorted range. Sort before calling.");
 
      checkIteratorsNotNULL(begin, end);
 
      const Size n = std::distance(begin, end);
      if (n == 0)
      {
        throw Exception::InvalidRange(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION);
      }
      if (q < 0.0 || q > 1.0)
      {
        throw Exception::InvalidValue(__FILE__, __LINE__, OPENMS_PRETTY_FUNCTION,
                                      "q must be in [0,1]",StringUtils::toStr(q));
      }
      if (n == 1) return static_cast<double>(*begin);
 
      const double pos = q * static_cast<double>(n - 1);
      const Size i = static_cast<Size>(std::floor(pos));
      const double frac = pos - static_cast<double>(i);
 
      const auto it_i = begin + static_cast<typename std::iterator_traits<IteratorType>::difference_type>(i);
      if (frac == 0.0) return static_cast<double>(*it_i);
 
      const auto it_ip1 = it_i + 1;
      return (1.0 - frac) * static_cast<double>(*it_i) + frac * static_cast<double>(*it_ip1);
    }
 
    template <typename IteratorType>
    double tukeyUpperFence(IteratorType begin, IteratorType end, double k = 1.5)
    {
        std::vector<double> v;
        v.reserve(std::distance(begin, end));
        for (auto it = begin; it != end; ++it)
        {
          if (std::isfinite(*it)) v.push_back(static_cast<double>(*it));
        }
        if (v.size() < 4) return std::numeric_limits<double>::infinity();
 
        std::sort(v.begin(), v.end());
        const double q1  = quantile(v.begin(), v.end(), 0.25);
        const double q3  = quantile(v.begin(), v.end(), 0.75);
        const double iqr = q3 - q1;
        if (!(iqr > 0.0)) return std::numeric_limits<double>::infinity();
 
        return q3 + k * iqr;
    }
 
    template <typename IteratorType>
    double tailFractionAbove(IteratorType begin, IteratorType end, double threshold)
    {
        size_t n = 0, n_tail = 0;
        for (auto it = begin; it != end; ++it)
        {
          const double x = static_cast<double>(*it);
          if (!std::isfinite(x)) continue;
          ++n;
          if (x > threshold) ++n_tail;
        }
        return (n == 0) ? 0.0 : static_cast<double>(n_tail) / static_cast<double>(n);
    }
 
    template <typename IteratorType>
    double winsorizedQuantile(IteratorType begin, IteratorType end, double q, double upper_fence)
    {
        std::vector<double> v;
        v.reserve(std::distance(begin, end));
        for (auto it = begin; it != end; ++it)
        {
          const double x = static_cast<double>(*it);
          if (!std::isfinite(x)) continue;
          v.push_back(x);
        }
        if (v.empty()) return 0.0;
 
        if (std::isfinite(upper_fence))
        {
          for (double& x : v)
          {
            if (x > upper_fence) x = upper_fence;
            if (x < 0.0) x = 0.0; // defensive; useful when passing |residual|
          }
        }
        std::sort(v.begin(), v.end());
        return quantile(v.begin(), v.end(), q);
    }
 
    template <typename IteratorType>
    AdaptiveQuantileResult adaptiveQuantile(IteratorType begin, IteratorType end, double q,
                            double k = 1.5,
                            double r_sparse = 0.01,
                            double r_dense  = 0.10)
    {
        AdaptiveQuantileResult res;
 
        // Copy finite values
        std::vector<double> v;
        v.reserve(std::distance(begin, end));
        for (auto it = begin; it != end; ++it)
        {
          if (std::isfinite(*it)) v.push_back(static_cast<double>(*it));
        }
        if (v.empty())
        {
          return res;
        }
 
        std::sort(v.begin(), v.end());
        const double half_raw = quantile(v.begin(), v.end(), q);
 
        // Robust path (winsorization at Tukey fence)
        const double uf       = tukeyUpperFence(v.begin(), v.end(), k);
        const double r        = std::isfinite(uf) ? tailFractionAbove(v.begin(), v.end(), uf) : 0.0;
        const double half_rob = winsorizedQuantile(v.begin(), v.end(), q, uf);
 
        // Blend weight w(r)
        double w = 0.0;
        if (r_dense <= r_sparse)
        {
          w = (r > r_sparse) ? 1.0 : 0.0;
        }
        else
        {
          const double t = (r - r_sparse) / (r_dense - r_sparse);
          w = std::max(0.0, std::min(1.0, t));
        }
 
        res.half_raw      = half_raw;
        res.half_rob      = half_rob;
        res.upper_fence   = uf;
        res.tail_fraction = r;
        res.weight        = w;
        res.blended       = (1.0 - w) * half_rob + w * half_raw;
        return res;
    }
 
    template <typename IteratorType>
    static double variance(IteratorType begin, IteratorType end,
                           double mean = std::numeric_limits<double>::max())
    {
      checkIteratorsNotNULL(begin, end);
      double sum_value = 0.0;
      if (mean == std::numeric_limits<double>::max())
      {
        mean = Math::mean(begin, end);
      }
      for (IteratorType iter=begin; iter!=end; ++iter)
      {
        double diff = *iter - mean;
        sum_value += diff * diff;
      }
      return sum_value / (std::distance(begin, end)-1);
    }
 
    template <typename IteratorType>
    static double sd(IteratorType begin, IteratorType end,
                     double mean = std::numeric_limits<double>::max())
    {
      checkIteratorsNotNULL(begin, end);
      return std::sqrt( variance(begin, end, mean) );
    }
 
    template <typename IteratorType>
    static double absdev(IteratorType begin, IteratorType end,
                         double mean = std::numeric_limits<double>::max())
    {
      checkIteratorsNotNULL(begin, end);
      if (mean == std::numeric_limits<double>::max())
      {
        mean = Math::mean(begin, end);
      }
      return MeanAbsoluteDeviation(begin, end, mean);
    }
 
    template <typename IteratorType1, typename IteratorType2>
    static double covariance(IteratorType1 begin_a, IteratorType1 end_a,
                             IteratorType2 begin_b, IteratorType2 end_b)
    {
      //no data or different lengths
      checkIteratorsNotNULL(begin_a, end_a);
 
      double sum_value = 0.0;
      double mean_a = Math::mean(begin_a, end_a);
      double mean_b = Math::mean(begin_b, end_b);
      IteratorType1 iter_a = begin_a;
      IteratorType2 iter_b = begin_b;
      for (; iter_a != end_a; ++iter_a, ++iter_b)
      {
        /* assure both ranges have the same number of elements */
        checkIteratorsAreValid(begin_b, end_b, begin_a, end_a);
        sum_value += (*iter_a - mean_a) * (*iter_b - mean_b);
      }
      /* assure both ranges have the same number of elements */
      checkIteratorsEqual(iter_b, end_b);
      Size n = std::distance(begin_a, end_a);
      return sum_value / (n-1);
    }
 
    template <typename IteratorType1, typename IteratorType2>
    static double meanSquareError(IteratorType1 begin_a, IteratorType1 end_a,
                                  IteratorType2 begin_b, IteratorType2 end_b)
    {
      //no data or different lengths
      checkIteratorsNotNULL(begin_a, end_a);
 
      SignedSize dist = std::distance(begin_a, end_a);
      double error = 0;
      IteratorType1 iter_a = begin_a;
      IteratorType2 iter_b = begin_b;
      for (; iter_a != end_a; ++iter_a, ++iter_b)
      {
        /* assure both ranges have the same number of elements */
        checkIteratorsAreValid(iter_b, end_b, iter_a, end_a);
 
        double tmp(*iter_a - *iter_b);
        error += tmp * tmp;
      }
      /* assure both ranges have the same number of elements */
      checkIteratorsEqual(iter_b, end_b);
 
      return error / dist;
    }
 
    template <typename IteratorType1, typename IteratorType2>
    static double rootMeanSquareError(IteratorType1 begin_a, IteratorType1 end_a,
                                      IteratorType2 begin_b, IteratorType2 end_b)
    {
      return std::sqrt(meanSquareError(begin_a, end_a, begin_b, end_b));
    }
 
    template <typename IteratorType1, typename IteratorType2>
    static double classificationRate(IteratorType1 begin_a, IteratorType1 end_a,
                                     IteratorType2 begin_b, IteratorType2 end_b)
    {
      //no data or different lengths
      checkIteratorsNotNULL(begin_a, end_a);
 
      SignedSize dist = std::distance(begin_a, end_a);
      SignedSize correct = dist;
      IteratorType1 iter_a = begin_a;
      IteratorType2 iter_b = begin_b;
      for (; iter_a != end_a; ++iter_a, ++iter_b)
      {
        /* assure both ranges have the same number of elements */
        checkIteratorsAreValid(iter_b, end_b, iter_a, end_a);
        if ((*iter_a < 0 && *iter_b >= 0) || (*iter_a >= 0 && *iter_b < 0))
        {
          --correct;
        }
 
      }
      /* assure both ranges have the same number of elements */
      checkIteratorsEqual(iter_b, end_b);
 
      return double(correct) / dist;
    }
 
    template <typename IteratorType1, typename IteratorType2>
    static double matthewsCorrelationCoefficient(
      IteratorType1 begin_a, IteratorType1 end_a,
      IteratorType2 begin_b, IteratorType2 end_b)
    {
      //no data or different lengths
      checkIteratorsNotNULL(begin_a, end_b);
 
      double tp = 0;
      double fp = 0;
      double tn = 0;
      double fn = 0;
      IteratorType1 iter_a = begin_a;
      IteratorType2 iter_b = begin_b;
      for (; iter_a != end_a; ++iter_a, ++iter_b)
      {
        /* assure both ranges have the same number of elements */
        checkIteratorsAreValid(iter_b, end_b, iter_a, end_a);
 
        if (*iter_a < 0 && *iter_b >= 0)
        {
          ++fn;
        }
        else if (*iter_a < 0 && *iter_b < 0)
        {
          ++tn;
        }
        else if (*iter_a >= 0 && *iter_b >= 0)
        {
          ++tp;
        }
        else if (*iter_a >= 0 && *iter_b < 0)
        {
          ++fp;
        }
      }
      /* assure both ranges have the same number of elements */
      checkIteratorsEqual(iter_b, end_b);
 
      return (tp * tn - fp * fn) / std::sqrt((tp + fp) * (tp + fn) * (tn + fp) * (tn + fn));
    }
 
    template <typename IteratorType1, typename IteratorType2>
    static double pearsonCorrelationCoefficient(
      IteratorType1 begin_a, IteratorType1 end_a,
      IteratorType2 begin_b, IteratorType2 end_b)
    {
      //no data or different lengths
      checkIteratorsNotNULL(begin_a, end_a);
 
      //calculate average
      SignedSize dist = std::distance(begin_a, end_a);
      double avg_a = std::accumulate(begin_a, end_a, 0.0) / dist;
      double avg_b = std::accumulate(begin_b, end_b, 0.0) / dist;
 
      double numerator = 0;
      double denominator_a = 0;
      double denominator_b = 0;
      IteratorType1 iter_a = begin_a;
      IteratorType2 iter_b = begin_b;
      for (; iter_a != end_a; ++iter_a, ++iter_b)
      {
        /* assure both ranges have the same number of elements */
        checkIteratorsAreValid(iter_b, end_b, iter_a, end_a);
        double temp_a = *iter_a - avg_a;
        double temp_b = *iter_b - avg_b;
        numerator += (temp_a * temp_b);
        denominator_a += (temp_a * temp_a);
        denominator_b += (temp_b * temp_b);
      }
      /* assure both ranges have the same number of elements */
      checkIteratorsEqual(iter_b, end_b);
      return numerator / std::sqrt(denominator_a * denominator_b);
    }
 
    template <typename Value>
    static void computeRank(std::vector<Value> & w)
    {
      Size i = 0; // main index
      Size z  = 0;  // "secondary" index
      Value rank = 0;
      Size n = (w.size() - 1);
      //store original indices for later
      std::vector<std::pair<Size, Value> > w_idx;
      for (Size j = 0; j < w.size(); ++j)
      {
        w_idx.push_back(std::make_pair(j, w[j]));
      }
      //sort
      std::sort(w_idx.begin(), w_idx.end(),
                [](const auto& pair1, const auto& pair2) { return pair1.second < pair2.second; });
      //replace pairs <orig_index, value> in w_idx by pairs <orig_index, rank>
      while (i < n)
      {
        // test for equality with tolerance:
        if (fabs(w_idx[i + 1].second - w_idx[i].second) > 0.0000001 * fabs(w_idx[i + 1].second)) // no tie
        {
          w_idx[i].second = Value(i + 1);
          ++i;
        }
        else // tie, replace by mean rank
        {
          // count number of ties
          for (z = i + 1; (z <= n) && fabs(w_idx[z].second - w_idx[i].second) <= 0.0000001 * fabs(w_idx[z].second); ++z)
          {
          }
          // compute mean rank of tie
          rank = 0.5 * (i + z + 1);
          // replace intensities by rank
          for (Size v = i; v <= z - 1; ++v)
          {
            w_idx[v].second = rank;
          }
          i = z;
        }
      }
      if (i == n)
        w_idx[n].second = Value(n + 1);
      //restore original order and replace elements of w with their ranks
      for (Size j = 0; j < w.size(); ++j)
      {
        w[w_idx[j].first] = w_idx[j].second;
      }
    }
 
    template <typename IteratorType1, typename IteratorType2>
    static double rankCorrelationCoefficient(
      IteratorType1 begin_a, IteratorType1 end_a,
      IteratorType2 begin_b, IteratorType2 end_b)
    {
      //no data or different lengths
      checkIteratorsNotNULL(begin_a, end_a);
 
      // store and sort intensities of model and data
      SignedSize dist = std::distance(begin_a, end_a);
      std::vector<double> ranks_data;
      ranks_data.reserve(dist);
      std::vector<double> ranks_model;
      ranks_model.reserve(dist);
      IteratorType1 iter_a = begin_a;
      IteratorType2 iter_b = begin_b;
      for (; iter_a != end_a; ++iter_a, ++iter_b)
      {
        /* assure both ranges have the same number of elements */
        checkIteratorsAreValid(iter_b, end_b, iter_a, end_a);
 
        ranks_model.push_back(*iter_a);
        ranks_data.push_back(*iter_b);
      }
      /* assure both ranges have the same number of elements */
      checkIteratorsEqual(iter_b, end_b);
 
      // replace entries by their ranks
      computeRank(ranks_data);
      computeRank(ranks_model);
 
      double mu = double(ranks_data.size() + 1) / 2.; // mean of ranks
      // Was the following, but I think the above is more correct ... (Clemens)
      // double mu = (ranks_data.size() + 1) / 2;
 
      double sum_model_data = 0;
      double sqsum_data = 0;
      double sqsum_model = 0;
 
      for (Int i = 0; i < dist; ++i)
      {
        sum_model_data += (ranks_data[i] - mu) * (ranks_model[i] - mu);
        sqsum_data += (ranks_data[i] - mu) * (ranks_data[i] - mu);
        sqsum_model += (ranks_model[i] - mu) * (ranks_model[i] - mu);
      }
 
      // check for division by zero
      if (!sqsum_data || !sqsum_model)
      {
        return 0;
      }
 
      return sum_model_data / (std::sqrt(sqsum_data) * std::sqrt(sqsum_model));
    }
 
    template<typename T>
    struct SummaryStatistics
    {
      SummaryStatistics() = default;
 
      // Ctor with data
      SummaryStatistics(T& data)
      {
        count = data.size();
        // Sanity check: avoid core dump if no data points present.
        if (data.empty())
        {
          mean = variance = min = lowerq = median = upperq = max = 0.0;
        }
        else
        {
          sort(data.begin(), data.end());
          mean = Math::mean(data.begin(), data.end());
          // variance divides by (n-1) and is undefined for a single value; report 0.0 as in the empty case
          variance = (count > 1) ? Math::variance(data.begin(), data.end(), mean) : 0.0;
          min = data.front();
          lowerq = Math::quantile1st(data.begin(), data.end(), true);
          median = Math::median(data.begin(), data.end(), true);
          upperq = Math::quantile3rd(data.begin(), data.end(), true);
          max = data.back();
        }
      }
 
      double mean = 0, variance = 0 , lowerq = 0, median = 0, upperq = 0;
      typename T::value_type min = 0, max = 0;
      size_t count = 0;
    };
 
  }   // namespace Math
} // namespace OpenMS