source: branches/kendall-score/yat/statistics/Kendall.cc @ 4068

// $Id: Kendall.cc 4068 2021-08-05 23:28:50Z peter $

/*
  Copyright (C) 2011, 2012, 2020 Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 3 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "Kendall.h"

#include <yat/utility/Ranking.h>
#include <yat/utility/stl_utility.h>

#include <gsl/gsl_cdf.h>

#include <boost/scoped_ptr.hpp>

#include <algorithm>
#include <cassert>
#include <cmath>
#include <iterator>
#include <limits>
#include <map>
#include <set>
#include <utility>
#include <vector>


#include <iostream> // debug
namespace theplu {
namespace yat {
namespace statistics {


  /**
     Calculate sum over all pair

     \sum_ij f((first[i]-first[j])(first2[i]-first2[j]))
     where f(x) is -1 if x<0,
                    0 if x=0,
               and +1 if x>0.
   */
  template<typename Iterator1, typename Iterator2>
  long int count(Iterator1 first1, Iterator1 last1, Iterator2 first2)
  {
    long int count=0;
    for (Iterator1 i=first1 ; i!=last1; ++i)
      for (Iterator1 j=first1; j<i; ++j) {
        if (*i==*j || first2[i-first1]==first2[j-first1])
          continue;
        if ((*i > *j) == (first2[i-first1] > first2[j-first1]))
          ++count;
        else
          --count;
      }
    return count;
  }



  class Kendall::Pimpl
  {
    class Count
    {
    public:
      Count(const std::multiset<std::pair<double, double> >& data);
      long int count(void) const;
      double score(void) const;
      class Ties
      {
      public:
        Ties(void);
        void add(size_t n);
        bool have_ties(void) const { return n_pairs_; }
        // \return \sum x * (x-1)
        unsigned long int n_pairs(void) const { return n_pairs_;}
        // \return \sum x * (x-1) * (x-2)
        unsigned long int n_triples(void) const { return n_triples_; }
        // \return \sum x * (x-1) * (2*x+5)
        unsigned long int v_correction(void) const { return v_correction_; }
      private:
        unsigned long int n_pairs_;
        unsigned long int n_triples_;
        unsigned long int v_correction_;
      };

      const Ties& x_ties(void) const;
      const Ties& y_ties(void) const;

    private:
      Ties x_ties_;
      Ties y_ties_;

      // # pairs such that (x_i-x_j)(y_i-y_j) > 0
      long int concordant_;
      // # pairs such that (x_i-x_j)(y_i-y_j) < 0
      long int discordant_;
      // # pairs such that x_i!=x_j && y_i==y_j
      long int extraX_;
      // # pairs such that x_i==x_j && y_i!=y_j
      long int extraY_;
      // # pairs such that x_i==x_j && y_i==y_j
      //long int spare_;

      template<typename Iterator>
      void calculate_ties(Iterator first, Iterator last, Ties& ties)
      {
        while (first != last) {
          Iterator upper = first;
          size_t n = 1;
          ++upper;
          while (upper!=last && *upper==*first) {
            ++n;
            ++upper;
          }
          ties.add(n);
          first = upper;
        }
      }
    };

  public:
    Pimpl(void);
    Pimpl(const Pimpl& other);
    Pimpl& operator=(const Pimpl& rhs);
    void add(double x, double y);
    size_t n(void) const;
    /// \return one-sided p-value
    double p_approx(bool right) const;
    double p_exact(bool right, bool left) const;
    void reset(void);
    double score(void) const;
  private:
    // return # concordant pairs minus # discordant pairs
    long int count(void) const;
    // return estimated variance of score
    double variance(void) const;
    // data always sort wrt first and then second (if first are equal)
    std::multiset<std::pair<double, double> > data_;
    // calculated in score(void)
    boost::scoped_ptr<Count> count_;
  };


  // Kendall class
  Kendall::Kendall(void)
    : pimpl_(new Pimpl)
  {
  }


  Kendall::Kendall(const Kendall& rhs)
    : pimpl_(new Pimpl(*rhs.pimpl_))
  {
  }


  Kendall::Kendall(Kendall&& rhs)
    : pimpl_(rhs.pimpl_)
  {
    rhs.pimpl_ = nullptr;
  }


  Kendall::~Kendall(void)
  {
    delete pimpl_;
  }


  void Kendall::add(double x, double y)
  {
    pimpl_->add(x, y);
  }


  size_t Kendall::n(void) const
  {
    return pimpl_->n();
  }


  double Kendall::score(void) const
  {
    return pimpl_->score();
  }


  double Kendall::p_left(bool exact) const
  {
    if (!exact)
      return pimpl_->p_approx(false);
    return pimpl_->p_exact(false, true);
  }


  double Kendall::p_right(bool exact) const
  {
    if (!exact)
      return pimpl_->p_approx(true);
    return pimpl_->p_exact(true, false);
  }


  double Kendall::p_value(bool exact) const
  {
    if (exact)
      return pimpl_->p_exact(true, true);
    if (score()>0.0)
      return 2*p_right(false);
    return 2*p_left(false);
  }


  void Kendall::reset(void)
  {
    pimpl_->reset();
  }


  Kendall& Kendall::operator=(const Kendall& rhs)
  {
    if (&rhs == this)
      return *this;

    assert(pimpl_);
    assert(rhs.pimpl_);
    *pimpl_ = *rhs.pimpl_;
    return *this;
  }


  Kendall& Kendall::operator=(Kendall&& rhs)
  {
    std::swap(pimpl_, rhs.pimpl_);
    return *this;
  }


  Kendall::Pimpl::Count::Count(const std::multiset<std::pair<double,double>>& data)
  {
    // We follow 3 Algorithm SDTau for some-duplicate datasets in
    // 'Fast Algorithms For The Calculation Of Kendall's Tau'
    // by David Christen (Computational Statistics, March 2005)

    // data is sorted w.r.t. ::first
    calculate_ties(utility::pair_first_iterator(data.begin()),
                   utility::pair_first_iterator(data.end()),
                   x_ties_);

    /*
                y1 < y2  y2 == y2  y2 > y2
      x1 <  x2     C        eX        D
      x1 == x2    eY       spare      -
      x1 >  x2     -        -         -

      We categorise pairs into five categories:
      C: Concordant
      D: Discordant
      eX: extra X; Ys and only Ys are equal
      eY: extra Y; Xs and only Xs are equal
      spare: both Xs and Yy are equal

      Due to symmetry reasons and because data container is sorted, we
      can ignore lower part of the matrix above.
     */

    concordant_ = 0;
    discordant_ = 0;
    extraX_ = 0;
    extraY_ = 0;

    unsigned long int eY = 0;
    // size of the current equal range, i.e., number of data points
    // for X_i : X_j == X_i, Y_j == Y_i, j <= i including the current
    // point
    unsigned long int ties = 1; // because loop below skip first entry
    utility::Ranking<double> Y;

    // loop over data, which is sorted w.r.t. ::first
    auto previous = data.cbegin();
    assert(previous != data.cend());
    Y.insert(previous->second);
    auto it = std::next(previous);
    while (it!=data.cend()) {
      assert(previous->first <= it->first);
      // X not equal
      if (it->first != previous->first) {
        eY = 0;
        ties = 1;
      }
      // y also equal
      else if (it->second == previous->second)
        ++ties;
      else { // x equal, y not equal
        eY += ties;
        ties = 1;
      }

      Y.insert(it->second);
      // FIXME can we use return value from insert instead
      auto lower = Y.lower_bound(it->second);
      // number of element in Y smaller than it->second
      int n_smaller = Y.ranking(lower);
      // number of element in Y equal to it->second
      int n_equal = 1;
      assert(lower != Y.cend());
      auto upper = std::next(lower);
      while (upper!=Y.cend() && *upper==*lower) {
        ++upper;
        ++n_equal;
      }
      size_t i = Y.size();

      // n_smaller (y<yi) is the union of concordant (y<yi,x<xi)
      // and eY (y<yi,x==xi)
      int C = n_smaller - eY;

      int eX =  n_equal - ties;

      int D = i - (C + eX + eY + ties);

      extraY_ += eY;
      extraX_ += eX;
      concordant_ += C;
      discordant_ += D;
      previous = it;
      ++it;
    }

  }


  long int Kendall::Pimpl::Count::count(void) const
  {
    return concordant_ - discordant_;
  }


  double Kendall::Pimpl::Count::score(void) const
  {
    double numerator = count();
    double denominator = concordant_ + discordant_;
    if (extraX_ || extraY_) {
      denominator =
        std::sqrt((denominator + extraX_)*(denominator + extraY_));
    }
    return numerator / denominator;
  }


  const Kendall::Pimpl::Count::Ties& Kendall::Pimpl::Count::x_ties(void) const
  {
    return x_ties_;
  }


  const Kendall::Pimpl::Count::Ties& Kendall::Pimpl::Count::y_ties(void) const
  {
    return y_ties_;
  }


  double Kendall::Pimpl::variance(void) const
  {
    /*
      According to wikipedia,
      z = k / sqrt(v)
      is approximately standard normal
      v = (v0 - vt - vu)/18 + v1 + v2
      v0 = n(n-1)(2n+5)
      vt = \sum t(t-1)(2t+5)
      vu = \sum u(u-1)(2u+5)
      v1 = \sum t(t-1)) * \sum u(u-1) / (2n(n-1))
      v2 = sum t(t-1)(t-2) \sum u(u-1)(u-2) / (9n(n-1)(n-2))

      where t is number of equal values in group i and similarly u for
      y.
    */
    double n = data_.size();
    double v0 = n*(n-1)*(2*n+5);
    double vt = 0;
    double vu = 0;
    double v1 = 0;
    double v2 = 0;
    assert(count_);
    auto& x_ties = count_->x_ties();
    auto& y_ties = count_->y_ties();
    // all correction terms above are zero in absence of ties
    bool x_have_ties = x_ties.have_ties();
    bool y_have_ties = y_ties.have_ties();
    if (x_have_ties || y_have_ties) {
      if (x_have_ties)
        vt = x_ties.v_correction();
      if (y_have_ties) {
        vu = y_ties.v_correction();
        if (x_have_ties) {
          v1 = x_ties.n_pairs() * (y_ties.n_pairs() / (2*n*(n-1)));
          v2 = x_ties.n_triples();
          if (v2)
            v2 *= y_ties.n_triples() / (9*n*(n-1)*(n-2));
        }
      }
    }
    return (v0 - vt - vu)/18 + v1 + v2;
  }


  Kendall::Pimpl::Count::Ties::Ties(void)
    : n_pairs_(0), n_triples_(0), v_correction_(0)
  {}


  void Kendall::Pimpl::Count::Ties::add(size_t n)
  {
    unsigned long int factor = n * (n-1);
    n_pairs_ += factor;
    n_triples_ += factor * (n-2);
    v_correction_ += factor * (2*n+5);
  }


  Kendall::Pimpl::Pimpl(void)
  {}


  Kendall::Pimpl::Pimpl(const Pimpl& other)
    : data_(other.data_)
  {}


  Kendall::Pimpl& Kendall::Pimpl::operator=(const Pimpl& rhs)
  {
    data_ = rhs.data_;
    count_.reset();
    return *this;
  }


  void Kendall::Pimpl::add(double x, double y)
  {
    data_.insert(std::make_pair(x, y));
    count_.reset();
  }


  size_t Kendall::Pimpl::n(void) const
  {
    return data_.size();
  }


  double Kendall::Pimpl::p_approx(bool right) const
  {
    double k = count();
    if (!right)
      k = -k;
    return gsl_cdf_gaussian_Q(k, std::sqrt(variance()));
  }


  double Kendall::Pimpl::p_exact(bool right, bool left) const
  {
    long int upper = 0;
    long int lower = 0;
    if (right) {
      if (left) {
        upper = std::max(count(), -count());
        lower = -upper;
      }
      else {
        upper = count();
        lower = std::numeric_limits<long int>::min();
      }
    }
    else {
      assert(left && "left or right must be true");
      upper = std::numeric_limits<long int>::max();
      lower = count();
    }

    // create a copy of the data, sort it with respect to ::second and
    // then iterate through the permutations of second while keeping
    // first constant. It means we need to do one extra initial sort,
    // but OTOH the permuted data is always almost sorted.
    std::vector<std::pair<double,double>> data(data_.begin(), data_.end());
    using utility::pair_second_iterator;
    std::sort(pair_second_iterator(data.begin()),
              pair_second_iterator(data.end()));
    unsigned int n = 0;
    unsigned int total = 0;
    do {
      std::multiset<std::pair<double,double>>
                    dataset(data.begin(), data.end());
      Count count(dataset);
      if (count.count() <= lower || count.count() >= upper)
        ++n;
      ++total;
    }
    while (std::next_permutation(pair_second_iterator(data.begin()),
                                 pair_second_iterator(data.end())));

    return static_cast<double>(n)/static_cast<double>(total);
  }


  void Kendall::Pimpl::reset(void)
  {
    Pimpl tmp;
    *this = tmp;
  }


  double Kendall::Pimpl::score(void) const
  {
    count();
    assert(count_.get());
    return count_->score();
  }


  long int Kendall::Pimpl::count(void) const
  {
    if (!count_)
      // const_cast to allow lazy eval is more restrictive than
      // making count_ mutable.
      const_cast<Pimpl*>(this)->count_.reset(new Count(data_));
    return count_->count();
  }

}}} // of namespace statistics, yat, and theplu