source: trunk/yat/statistics/Kendall.cc @ 4117

Last change on this file since 4117 was 4117, checked in by Peter, 7 months ago

add assert

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 12.1 KB
Line 
1// $Id: Kendall.cc 4117 2021-10-22 12:24:14Z peter $
2
3/*
4  Copyright (C) 2011, 2012, 2020, 2021 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 <utility>
41#include <vector>
42
43namespace theplu {
44namespace yat {
45namespace statistics {
46
47
48  /**
49     Calculate sum over all pair
50
51     \sum_ij f((first[i]-first[j])(first2[i]-first2[j]))
52     where f(x) is -1 if x<0,
53                    0 if x=0,
54               and +1 if x>0.
55   */
56  template<typename Iterator1, typename Iterator2>
57  long int count(Iterator1 first1, Iterator1 last1, Iterator2 first2)
58  {
59    long int count=0;
60    for (Iterator1 i=first1 ; i!=last1; ++i)
61      for (Iterator1 j=first1; j<i; ++j) {
62        if (*i==*j || first2[i-first1]==first2[j-first1])
63          continue;
64        if ((*i > *j) == (first2[i-first1] > first2[j-first1]))
65          ++count;
66        else
67          --count;
68      }
69    return count;
70  }
71
72
73
74  class Kendall::Pimpl
75  {
76    class Count
77    {
78    public:
79      Count(const std::multiset<std::pair<double, double> >& data);
80      long int count(void) const;
81      double score(void) const;
82      class Ties
83      {
84      public:
85        Ties(void);
86        void add(size_t n);
87        bool have_ties(void) const { return n_pairs_; }
88        // \return \sum x * (x-1)
89        unsigned long int n_pairs(void) const { return n_pairs_;}
90        // \return \sum x * (x-1) * (x-2)
91        unsigned long int n_triples(void) const { return n_triples_; }
92        // \return \sum x * (x-1) * (2*x+5)
93        unsigned long int v_correction(void) const { return v_correction_; }
94      private:
95        unsigned long int n_pairs_;
96        unsigned long int n_triples_;
97        unsigned long int v_correction_;
98      };
99
100      const Ties& x_ties(void) const;
101      const Ties& y_ties(void) const;
102
103    private:
104      Ties x_ties_;
105      Ties y_ties_;
106
107      // # pairs such that (x_i-x_j)(y_i-y_j) > 0
108      long int concordant_;
109      // # pairs such that (x_i-x_j)(y_i-y_j) < 0
110      long int discordant_;
111      // # pairs such that x_i!=x_j && y_i==y_j
112      long int extraX_;
113      // # pairs such that x_i==x_j && y_i!=y_j
114      long int extraY_;
115      // # pairs such that x_i==x_j && y_i==y_j
116      //long int spare_;
117
118      template<typename Iterator>
119      void calculate_ties(Iterator first, Iterator last, Ties& ties)
120      {
121        assert(std::is_sorted(first, last));
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(Kendall&& rhs)
173    : pimpl_(rhs.pimpl_)
174  {
175    rhs.pimpl_ = nullptr;
176  }
177
178
179  Kendall::~Kendall(void)
180  {
181    delete pimpl_;
182  }
183
184
185  void Kendall::add(double x, double y)
186  {
187    pimpl_->add(x, y);
188  }
189
190
191  size_t Kendall::n(void) const
192  {
193    return pimpl_->n();
194  }
195
196
197  double Kendall::score(void) const
198  {
199    return pimpl_->score();
200  }
201
202
203  double Kendall::p_left(bool exact) const
204  {
205    if (!exact)
206      return pimpl_->p_approx(false);
207    return pimpl_->p_exact(false, true);
208  }
209
210
211  double Kendall::p_right(bool exact) const
212  {
213    if (!exact)
214      return pimpl_->p_approx(true);
215    return pimpl_->p_exact(true, false);
216  }
217
218
219  double Kendall::p_value(bool exact) const
220  {
221    if (exact)
222      return pimpl_->p_exact(true, true);
223    if (score()>0.0)
224      return 2*p_right(false);
225    return 2*p_left(false);
226  }
227
228
229  void Kendall::reset(void)
230  {
231    pimpl_->reset();
232  }
233
234
235  Kendall& Kendall::operator=(const Kendall& rhs)
236  {
237    if (&rhs == this)
238      return *this;
239
240    assert(pimpl_);
241    assert(rhs.pimpl_);
242    *pimpl_ = *rhs.pimpl_;
243    return *this;
244  }
245
246
247  Kendall& Kendall::operator=(Kendall&& rhs)
248  {
249    std::swap(pimpl_, rhs.pimpl_);
250    return *this;
251  }
252
253
254  Kendall::Pimpl::Count::Count(const std::multiset<std::pair<double,double>>& data)
255  {
256    // We follow 3 Algorithm SDTau for some-duplicate datasets in
257    // 'Fast Algorithms For The Calculation Of Kendall's Tau'
258    // by David Christen (Computational Statistics, March 2005)
259
260    // data is sorted w.r.t. ::first
261    calculate_ties(utility::pair_first_iterator(data.begin()),
262                   utility::pair_first_iterator(data.end()),
263                   x_ties_);
264
265    std::vector<double> vec(utility::pair_second_iterator(data.begin()),
266                            utility::pair_second_iterator(data.end()));
267    std::sort(vec.begin(), vec.end());
268    calculate_ties(vec.begin(), vec.end(), y_ties_);
269
270    /*
271                y1 < y2  y2 == y2  y2 > y2
272      x1 <  x2     C        eX        D
273      x1 == x2    eY       spare      -
274      x1 >  x2     -        -         -
275
276      We categorise pairs into five categories:
277      C: Concordant
278      D: Discordant
279      eX: extra X; Ys and only Ys are equal
280      eY: extra Y; Xs and only Xs are equal
281      spare: both Xs and Yy are equal
282
283      Due to symmetry reasons and because data container is sorted, we
284      can ignore lower part of the matrix above.
285     */
286
287    concordant_ = 0;
288    discordant_ = 0;
289    extraX_ = 0;
290    extraY_ = 0;
291
292    unsigned long int eY = 0;
293    // size of the current equal range, i.e., number of data points
294    // for X_i : X_j == X_i, Y_j == Y_i, j <= i including the current
295    // point
296    unsigned long int ties = 1; // because loop below skip first entry
297    utility::Ranking<double> Y;
298
299    // loop over data, which is sorted w.r.t. ::first
300    auto previous = data.cbegin();
301    assert(previous != data.cend());
302    Y.insert(previous->second);
303    auto it = std::next(previous);
304    while (it!=data.cend()) {
305      assert(previous->first <= it->first);
306      // X not equal
307      if (it->first != previous->first) {
308        eY = 0;
309        ties = 1;
310      }
311      // y also equal
312      else if (it->second == previous->second)
313        ++ties;
314      else { // x equal, y not equal
315        eY += ties;
316        ties = 1;
317      }
318
319      auto lower = Y.insert(it->second);
320      if (lower!=Y.begin() && *std::prev(lower) == it->second)
321        lower = Y.lower_bound(it->second);
322      // number of element in Y smaller than it->second
323      int n_smaller = Y.ranking(lower);
324      // number of element in Y equal to it->second
325      int n_equal = 1;
326      assert(lower != Y.cend());
327      auto upper = std::next(lower);
328      while (upper!=Y.cend() && *upper==*lower) {
329        ++upper;
330        ++n_equal;
331      }
332      size_t i = Y.size();
333
334      // n_smaller (y<yi) is the union of concordant (y<yi,x<xi)
335      // and eY (y<yi,x==xi)
336      int C = n_smaller - eY;
337
338      int eX =  n_equal - ties;
339
340      int D = i - (C + eX + eY + ties);
341
342      extraY_ += eY;
343      extraX_ += eX;
344      concordant_ += C;
345      discordant_ += D;
346      previous = it;
347      ++it;
348    }
349
350  }
351
352
353  long int Kendall::Pimpl::Count::count(void) const
354  {
355    return concordant_ - discordant_;
356  }
357
358
359  double Kendall::Pimpl::Count::score(void) const
360  {
361    double numerator = count();
362    double denominator = concordant_ + discordant_;
363    if (extraX_ || extraY_) {
364      denominator =
365        std::sqrt((denominator + extraX_)*(denominator + extraY_));
366    }
367    return numerator / denominator;
368  }
369
370
371  const Kendall::Pimpl::Count::Ties& Kendall::Pimpl::Count::x_ties(void) const
372  {
373    return x_ties_;
374  }
375
376
377  const Kendall::Pimpl::Count::Ties& Kendall::Pimpl::Count::y_ties(void) const
378  {
379    return y_ties_;
380  }
381
382
383  double Kendall::Pimpl::variance(void) const
384  {
385    /*
386      According to wikipedia,
387      z = k / sqrt(v)
388      is approximately standard normal
389      v = (v0 - vt - vu)/18 + v1 + v2
390      v0 = n(n-1)(2n+5)
391      vt = \sum t(t-1)(2t+5)
392      vu = \sum u(u-1)(2u+5)
393      v1 = \sum t(t-1)) * \sum u(u-1) / (2n(n-1))
394      v2 = sum t(t-1)(t-2) \sum u(u-1)(u-2) / (9n(n-1)(n-2))
395
396      where t is number of equal x values in group i and similarly u for
397      y.
398    */
399    double n = data_.size();
400    double v0 = n*(n-1)*(2*n+5);
401    double vt = 0;
402    double vu = 0;
403    double v1 = 0;
404    double v2 = 0;
405    assert(count_);
406    auto& x_ties = count_->x_ties();
407    auto& y_ties = count_->y_ties();
408    // all correction terms above are zero in absence of ties
409    bool x_have_ties = x_ties.have_ties();
410    bool y_have_ties = y_ties.have_ties();
411    if (x_have_ties || y_have_ties) {
412      if (x_have_ties)
413        vt = x_ties.v_correction();
414      if (y_have_ties) {
415        vu = y_ties.v_correction();
416        if (x_have_ties) {
417          v1 = x_ties.n_pairs() * (y_ties.n_pairs() / (2*n*(n-1)));
418          v2 = x_ties.n_triples();
419          if (v2)
420            v2 *= y_ties.n_triples() / (9*n*(n-1)*(n-2));
421        }
422      }
423    }
424    return (v0 - vt - vu)/18 + v1 + v2;
425  }
426
427
428  Kendall::Pimpl::Count::Ties::Ties(void)
429    : n_pairs_(0), n_triples_(0), v_correction_(0)
430  {}
431
432
433  void Kendall::Pimpl::Count::Ties::add(size_t n)
434  {
435    unsigned long int factor = n * (n-1);
436    n_pairs_ += factor;
437    n_triples_ += factor * (n-2);
438    v_correction_ += factor * (2*n+5);
439  }
440
441
442  Kendall::Pimpl::Pimpl(void)
443  {}
444
445
446  Kendall::Pimpl::Pimpl(const Pimpl& other)
447    : data_(other.data_)
448  {}
449
450
451  Kendall::Pimpl& Kendall::Pimpl::operator=(const Pimpl& rhs)
452  {
453    data_ = rhs.data_;
454    count_.reset();
455    return *this;
456  }
457
458
459  void Kendall::Pimpl::add(double x, double y)
460  {
461    data_.insert(std::make_pair(x, y));
462    count_.reset();
463  }
464
465
466  size_t Kendall::Pimpl::n(void) const
467  {
468    return data_.size();
469  }
470
471
472  double Kendall::Pimpl::p_approx(bool right) const
473  {
474    double k = count();
475    if (!right)
476      k = -k;
477    return gsl_cdf_gaussian_Q(k, std::sqrt(variance()));
478  }
479
480
481  double Kendall::Pimpl::p_exact(bool right, bool left) const
482  {
483    long int upper = 0;
484    long int lower = 0;
485    if (right) {
486      if (left) {
487        upper = std::max(count(), -count());
488        lower = -upper;
489      }
490      else {
491        upper = count();
492        lower = std::numeric_limits<long int>::min();
493      }
494    }
495    else {
496      assert(left && "left or right must be true");
497      upper = std::numeric_limits<long int>::max();
498      lower = count();
499    }
500
501    // create a copy of the data, sort it with respect to ::second and
502    // then iterate through the permutations of second while keeping
503    // first constant. It means we need to do one extra initial sort,
504    // but OTOH the permuted data is always almost sorted.
505    std::vector<std::pair<double,double>> data(data_.begin(), data_.end());
506    using utility::pair_second_iterator;
507    std::sort(pair_second_iterator(data.begin()),
508              pair_second_iterator(data.end()));
509    unsigned int n = 0;
510    unsigned int total = 0;
511    do {
512      std::multiset<std::pair<double,double>>
513                    dataset(data.begin(), data.end());
514      Count count(dataset);
515      if (count.count() <= lower || count.count() >= upper)
516        ++n;
517      ++total;
518    }
519    while (std::next_permutation(pair_second_iterator(data.begin()),
520                                 pair_second_iterator(data.end())));
521
522    return static_cast<double>(n)/static_cast<double>(total);
523  }
524
525
526  void Kendall::Pimpl::reset(void)
527  {
528    Pimpl tmp;
529    *this = tmp;
530  }
531
532
533  double Kendall::Pimpl::score(void) const
534  {
535    count();
536    assert(count_.get());
537    return count_->score();
538  }
539
540
541  long int Kendall::Pimpl::count(void) const
542  {
543    if (!count_)
544      // const_cast to allow lazy eval is more restrictive than
545      // making count_ mutable.
546      const_cast<Pimpl*>(this)->count_.reset(new Count(data_));
547    return count_->count();
548  }
549
550}}} // of namespace statistics, yat, and theplu
Note: See TracBrowser for help on using the repository browser.