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

Last change on this file since 4064 was 4064, checked in by Peter, 12 months ago

first functioning version (tests pass) of Kendall class using the Ranking class. Split out the code in separate header and source files. Ranking class is still slow (linear) so Kendall::score is still quadratic. refs #710

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 11.6 KB
Line 
1// $Id: Kendall.cc 4064 2021-08-05 05:36:32Z peter $
2
3/*
4  Copyright (C) 2011, 2012, 2020 Peter Johansson
5
6  This file is part of the yat library, http://dev.thep.lu.se/yat
7
8  The yat library is free software; you can redistribute it and/or
9  modify it under the terms of the GNU General Public License as
10  published by the Free Software Foundation; either version 3 of the
11  License, or (at your option) any later version.
12
13  The yat library is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16  General Public License for more details.
17
18  You should have received a copy of the GNU General Public License
19  along with yat. If not, see <http://www.gnu.org/licenses/>.
20*/
21
22#include <config.h>
23
24#include "Kendall.h"
25
26#include <yat/utility/Ranking.h>
27#include <yat/utility/stl_utility.h>
28
29#include <gsl/gsl_cdf.h>
30
31#include <boost/scoped_ptr.hpp>
32
33#include <algorithm>
34#include <cassert>
35#include <cmath>
36#include <iterator>
37#include <limits>
38#include <map>
39#include <set>
40#include <vector>
41
42
43#include <iostream> // debug
44namespace theplu {
45namespace yat {
46namespace statistics {
47
48
49  /**
50     Calculate sum over all pair
51
52     \sum_ij f((first[i]-first[j])(first2[i]-first2[j]))
53     where f(x) is -1 if x<0,
54                    0 if x=0,
55               and +1 if x>0.
56   */
57  template<typename Iterator1, typename Iterator2>
58  long int count(Iterator1 first1, Iterator1 last1, Iterator2 first2)
59  {
60    long int count=0;
61    for (Iterator1 i=first1 ; i!=last1; ++i)
62      for (Iterator1 j=first1; j<i; ++j) {
63        if (*i==*j || first2[i-first1]==first2[j-first1])
64          continue;
65        if ((*i > *j) == (first2[i-first1] > first2[j-first1]))
66          ++count;
67        else
68          --count;
69      }
70    return count;
71  }
72
73
74
75  class Kendall::Pimpl
76  {
77    class Count
78    {
79    public:
80      Count(const std::multiset<std::pair<double, double> >& data);
81      long int count(void) const;
82      double score(void) const;
83      class Ties
84      {
85      public:
86        Ties(void);
87        void add(size_t n);
88        bool have_ties(void) const { return n_pairs_; }
89        // \return \sum x * (x-1)
90        unsigned long int n_pairs(void) const { return n_pairs_;}
91        // \return \sum x * (x-1) * (x-2)
92        unsigned long int n_triples(void) const { return n_triples_; }
93        // \return \sum x * (x-1) * (2*x+5)
94        unsigned long int v_correction(void) const { return v_correction_; }
95      private:
96        unsigned long int n_pairs_;
97        unsigned long int n_triples_;
98        unsigned long int v_correction_;
99      };
100
101      const Ties& x_ties(void) const;
102      const Ties& y_ties(void) const;
103
104    private:
105      Ties x_ties_;
106      Ties y_ties_;
107
108      // # pairs such that (x_i-x_j)(y_i-y_j) > 0
109      long int concordant_;
110      // # pairs such that (x_i-x_j)(y_i-y_j) < 0
111      long int discordant_;
112      // # pairs such that x_i!=x_j && y_i==y_j
113      long int extraX_;
114      // # pairs such that x_i==x_j && y_i!=y_j
115      long int extraY_;
116      // # pairs such that x_i==x_j && y_i==y_j
117      //long int spare_;
118
119      template<typename Iterator>
120      void calculate_ties(Iterator first, Iterator last, Ties& ties)
121      {
122        while (first != last) {
123          Iterator upper = first;
124          size_t n = 1;
125          ++upper;
126          while (upper!=last && *upper==*first) {
127            ++n;
128            ++upper;
129          }
130          ties.add(n);
131          first = upper;
132        }
133      }
134    };
135
136  public:
137    Pimpl(void);
138    Pimpl(const Pimpl& other);
139    Pimpl& operator=(const Pimpl& rhs);
140    void add(double x, double y);
141    size_t n(void) const;
142    /// \return one-sided p-value
143    double p_approx(bool right) const;
144    double p_exact(bool right, bool left) const;
145    void reset(void);
146    double score(void) const;
147  private:
148    // return # concordant pairs minus # discordant pairs
149    long int count(void) const;
150    // return estimated variance of score
151    double variance(void) const;
152    // data always sort wrt first and then second (if first are equal)
153    std::multiset<std::pair<double, double> > data_;
154    // calculated in score(void)
155    boost::scoped_ptr<Count> count_;
156  };
157
158
159  // Kendall class
160  Kendall::Kendall(void)
161    : pimpl_(new Pimpl)
162  {
163  }
164
165
166  Kendall::Kendall(const Kendall& rhs)
167    : pimpl_(new Pimpl(*rhs.pimpl_))
168  {
169  }
170
171
172  Kendall::~Kendall(void)
173  {
174    delete pimpl_;
175  }
176
177
178  void Kendall::add(double x, double y)
179  {
180    pimpl_->add(x, y);
181  }
182
183
184  size_t Kendall::n(void) const
185  {
186    return pimpl_->n();
187  }
188
189
190  double Kendall::score(void) const
191  {
192    return pimpl_->score();
193  }
194
195
196  double Kendall::p_left(bool exact) const
197  {
198    if (!exact)
199      return pimpl_->p_approx(false);
200    return pimpl_->p_exact(false, true);
201  }
202
203
204  double Kendall::p_right(bool exact) const
205  {
206    if (!exact)
207      return pimpl_->p_approx(true);
208    return pimpl_->p_exact(true, false);
209  }
210
211
212  double Kendall::p_value(bool exact) const
213  {
214    if (exact)
215      return pimpl_->p_exact(true, true);
216    if (score()>0.0)
217      return 2*p_right(false);
218    return 2*p_left(false);
219  }
220
221
222  void Kendall::reset(void)
223  {
224    pimpl_->reset();
225  }
226
227
228  Kendall& Kendall::operator=(const Kendall& rhs)
229  {
230    if (&rhs == this)
231      return *this;
232
233    assert(pimpl_);
234    assert(rhs.pimpl_);
235    *pimpl_ = *rhs.pimpl_;
236    return *this;
237  }
238
239
240  Kendall::Pimpl::Count::Count(const std::multiset<std::pair<double,double>>& data)
241  {
242    // We follow 3 Algorithm SDTau for some-duplicate datasets in
243    // 'Fast Algorithms For The Calculation Of Kendall's Tau'
244    // by David Christen (Computational Statistics, March 2005)
245
246    // data is sorted w.r.t. ::first
247    calculate_ties(utility::pair_first_iterator(data.begin()),
248                   utility::pair_first_iterator(data.end()),
249                   x_ties_);
250
251    /*
252                y1 < y2  y2 == y2  y2 > y2
253      x1 <  x2     C        eX        D
254      x1 == x2    eY       spare      -
255      x1 >  x2     -        -         -
256
257      We categorise pairs into five categories:
258      C: Concordant
259      D: Discordant
260      eX: extra X; Ys and only Ys are equal
261      eY: extra Y; Xs and only Xs are equal
262      spare: both Xs and Yy are equal
263
264      Due to symmetry reasons and because data container is sorted, we
265      can ignore lower part of the matrix above.
266     */
267
268    concordant_ = 0;
269    discordant_ = 0;
270    extraX_ = 0;
271    extraY_ = 0;
272
273    unsigned long int eY = 0;
274    // size of the current equal range, i.e., number of data points
275    // for X_i : X_j == X_i, Y_j == Y_i, j <= i including the current
276    // point
277    unsigned long int ties = 1; // because loop below skip first entry
278    utility::Ranking<double> Y;
279
280    // loop over data, which is sorted w.r.t. ::first
281    auto previous = data.cbegin();
282    assert(previous != data.cend());
283    Y.insert(previous->second);
284    auto it = std::next(previous);
285    while (it!=data.cend()) {
286      assert(previous->first <= it->first);
287      // X not equal
288      if (it->first != previous->first) {
289        eY = 0;
290        ties = 1;
291      }
292      // y also equal
293      else if (it->second == previous->second)
294        ++ties;
295      else { // x equal, y not equal
296        eY += ties;
297        ties = 1;
298      }
299
300      Y.insert(it->second);
301      // FIXME can we use return value from insert instead
302      auto lower = Y.lower_bound(it->second);
303      // number of element in Y smaller than it->second
304      int n_smaller = Y.ranking(lower);
305      // number of element in Y equal to it->second
306      int n_equal = 1;
307      assert(lower != Y.cend());
308      auto upper = std::next(lower);
309      while (upper!=Y.cend() && *upper==*lower) {
310        ++upper;
311        ++n_equal;
312      }
313      size_t i = Y.size();
314
315      // n_smaller (y<yi) is the union of concordant (y<yi,x<xi)
316      // and eY (y<yi,x==xi)
317      int C = n_smaller - eY;
318
319      int eX =  n_equal - ties;
320
321      int D = i - (C + eX + eY + ties);
322
323      extraY_ += eY;
324      extraX_ += eX;
325      concordant_ += C;
326      discordant_ += D;
327      previous = it;
328      ++it;
329    }
330
331  }
332
333
334  long int Kendall::Pimpl::Count::count(void) const
335  {
336    return concordant_ - discordant_;
337  }
338
339
340  double Kendall::Pimpl::Count::score(void) const
341  {
342    double numerator = count();
343    double denominator = concordant_ + discordant_;
344    if (extraX_ || extraY_) {
345      denominator =
346        std::sqrt((denominator + extraX_)*(denominator + extraY_));
347    }
348    return numerator / denominator;
349  }
350
351
352  const Kendall::Pimpl::Count::Ties& Kendall::Pimpl::Count::x_ties(void) const
353  {
354    return x_ties_;
355  }
356
357
358  const Kendall::Pimpl::Count::Ties& Kendall::Pimpl::Count::y_ties(void) const
359  {
360    return y_ties_;
361  }
362
363
364  double Kendall::Pimpl::variance(void) const
365  {
366    /*
367      According to wikipedia,
368      z = k / sqrt(v)
369      is approximately standard normal
370      v = (v0 - vt - vu)/18 + v1 + v2
371      v0 = n(n-1)(2n+5)
372      vt = \sum t(t-1)(2t+5)
373      vu = \sum u(u-1)(2u+5)
374      v1 = \sum t(t-1)) * \sum u(u-1) / (2n(n-1))
375      v2 = sum t(t-1)(t-2) \sum u(u-1)(u-2) / (9n(n-1)(n-2))
376
377      where t is number of equal values in group i and similarly u for
378      y.
379    */
380    double n = data_.size();
381    double v0 = n*(n-1)*(2*n+5);
382    double vt = 0;
383    double vu = 0;
384    double v1 = 0;
385    double v2 = 0;
386    assert(count_);
387    auto& x_ties = count_->x_ties();
388    auto& y_ties = count_->y_ties();
389    // all correction terms above are zero in absence of ties
390    bool x_have_ties = x_ties.have_ties();
391    bool y_have_ties = y_ties.have_ties();
392    if (x_have_ties || y_have_ties) {
393      if (x_have_ties)
394        vt = x_ties.v_correction();
395      if (y_have_ties) {
396        vu = y_ties.v_correction();
397        if (x_have_ties) {
398          v1 = x_ties.n_pairs() * (y_ties.n_pairs() / (2*n*(n-1)));
399          v2 = x_ties.n_triples();
400          if (v2)
401            v2 *= y_ties.n_triples() / (9*n*(n-1)*(n-2));
402        }
403      }
404    }
405    return (v0 - vt - vu)/18 + v1 + v2;
406  }
407
408
409  Kendall::Pimpl::Count::Ties::Ties(void)
410    : n_pairs_(0), n_triples_(0), v_correction_(0)
411  {}
412
413
414  void Kendall::Pimpl::Count::Ties::add(size_t n)
415  {
416    unsigned long int factor = n * (n-1);
417    n_pairs_ += factor;
418    n_triples_ += factor * (n-2);
419    v_correction_ += factor * (2*n+5);
420  }
421
422
423  Kendall::Pimpl::Pimpl(void)
424  {}
425
426
427  Kendall::Pimpl::Pimpl(const Pimpl& other)
428    : data_(other.data_)
429  {}
430
431
432  Kendall::Pimpl& Kendall::Pimpl::operator=(const Pimpl& rhs)
433  {
434    data_ = rhs.data_;
435    count_.reset();
436    return *this;
437  }
438
439
440  void Kendall::Pimpl::add(double x, double y)
441  {
442    data_.insert(std::make_pair(x, y));
443    count_.reset();
444  }
445
446
447  size_t Kendall::Pimpl::n(void) const
448  {
449    return data_.size();
450  }
451
452
453  double Kendall::Pimpl::p_approx(bool right) const
454  {
455    double k = count();
456    if (!right)
457      k = -k;
458    return gsl_cdf_gaussian_Q(k, std::sqrt(variance()));
459  }
460
461
462  double Kendall::Pimpl::p_exact(bool right, bool left) const
463  {
464    long int upper = 0;
465    long int lower = 0;
466    if (right) {
467      if (left) {
468        upper = std::max(count(), -count());
469        lower = -upper;
470      }
471      else {
472        upper = count();
473        lower = std::numeric_limits<long int>::min();
474      }
475    }
476    else {
477      assert(left && "left or right must be true");
478      upper = std::numeric_limits<long int>::max();
479      lower = count();
480    }
481
482    // create a copy of the data, sort it with respect to ::second and
483    // then iterate through the permutations of second while keeping
484    // first constant. It means we need to do one extra initial sort,
485    // but OTOH the permuted data is always almost sorted.
486    std::vector<std::pair<double,double>> data(data_.begin(), data_.end());
487    using utility::pair_second_iterator;
488    std::sort(pair_second_iterator(data.begin()),
489              pair_second_iterator(data.end()));
490    unsigned int n = 0;
491    unsigned int total = 0;
492    do {
493      std::multiset<std::pair<double,double>>
494                    dataset(data.begin(), data.end());
495      Count count(dataset);
496      if (count.count() <= lower || count.count() >= upper)
497        ++n;
498      ++total;
499    }
500    while (std::next_permutation(pair_second_iterator(data.begin()),
501                                 pair_second_iterator(data.end())));
502
503    return static_cast<double>(n)/static_cast<double>(total);
504  }
505
506
507  void Kendall::Pimpl::reset(void)
508  {
509    Pimpl tmp;
510    *this = tmp;
511  }
512
513
514  double Kendall::Pimpl::score(void) const
515  {
516    count();
517    assert(count_.get());
518    return count_->score();
519  }
520
521
522  long int Kendall::Pimpl::count(void) const
523  {
524    if (!count_)
525      // const_cast to allow lazy eval is more restrictive than
526      // making count_ mutable.
527      const_cast<Pimpl*>(this)->count_.reset(new Count(data_));
528    return count_->count();
529  }
530
531}}} // of namespace statistics, yat, and theplu
Note: See TracBrowser for help on using the repository browser.