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

Last change on this file since 4110 was 4110, checked in by Peter, 8 months ago

merge release 0.19 into trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 12.1 KB
Line 
1// $Id: Kendall.cc 4110 2021-09-27 23:35:15Z 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        while (first != last) {
122          Iterator upper = first;
123          size_t n = 1;
124          ++upper;
125          while (upper!=last && *upper==*first) {
126            ++n;
127            ++upper;
128          }
129          ties.add(n);
130          first = upper;
131        }
132      }
133    };
134
135  public:
136    Pimpl(void);
137    Pimpl(const Pimpl& other);
138    Pimpl& operator=(const Pimpl& rhs);
139    void add(double x, double y);
140    size_t n(void) const;
141    /// \return one-sided p-value
142    double p_approx(bool right) const;
143    double p_exact(bool right, bool left) const;
144    void reset(void);
145    double score(void) const;
146  private:
147    // return # concordant pairs minus # discordant pairs
148    long int count(void) const;
149    // return estimated variance of score
150    double variance(void) const;
151    // data always sort wrt first and then second (if first are equal)
152    std::multiset<std::pair<double, double> > data_;
153    // calculated in score(void)
154    boost::scoped_ptr<Count> count_;
155  };
156
157
158  // Kendall class
159  Kendall::Kendall(void)
160    : pimpl_(new Pimpl)
161  {
162  }
163
164
165  Kendall::Kendall(const Kendall& rhs)
166    : pimpl_(new Pimpl(*rhs.pimpl_))
167  {
168  }
169
170
171  Kendall::Kendall(Kendall&& rhs)
172    : pimpl_(rhs.pimpl_)
173  {
174    rhs.pimpl_ = nullptr;
175  }
176
177
178  Kendall::~Kendall(void)
179  {
180    delete pimpl_;
181  }
182
183
184  void Kendall::add(double x, double y)
185  {
186    pimpl_->add(x, y);
187  }
188
189
190  size_t Kendall::n(void) const
191  {
192    return pimpl_->n();
193  }
194
195
196  double Kendall::score(void) const
197  {
198    return pimpl_->score();
199  }
200
201
202  double Kendall::p_left(bool exact) const
203  {
204    if (!exact)
205      return pimpl_->p_approx(false);
206    return pimpl_->p_exact(false, true);
207  }
208
209
210  double Kendall::p_right(bool exact) const
211  {
212    if (!exact)
213      return pimpl_->p_approx(true);
214    return pimpl_->p_exact(true, false);
215  }
216
217
218  double Kendall::p_value(bool exact) const
219  {
220    if (exact)
221      return pimpl_->p_exact(true, true);
222    if (score()>0.0)
223      return 2*p_right(false);
224    return 2*p_left(false);
225  }
226
227
228  void Kendall::reset(void)
229  {
230    pimpl_->reset();
231  }
232
233
234  Kendall& Kendall::operator=(const Kendall& rhs)
235  {
236    if (&rhs == this)
237      return *this;
238
239    assert(pimpl_);
240    assert(rhs.pimpl_);
241    *pimpl_ = *rhs.pimpl_;
242    return *this;
243  }
244
245
246  Kendall& Kendall::operator=(Kendall&& rhs)
247  {
248    std::swap(pimpl_, rhs.pimpl_);
249    return *this;
250  }
251
252
253  Kendall::Pimpl::Count::Count(const std::multiset<std::pair<double,double>>& data)
254  {
255    // We follow 3 Algorithm SDTau for some-duplicate datasets in
256    // 'Fast Algorithms For The Calculation Of Kendall's Tau'
257    // by David Christen (Computational Statistics, March 2005)
258
259    // data is sorted w.r.t. ::first
260    calculate_ties(utility::pair_first_iterator(data.begin()),
261                   utility::pair_first_iterator(data.end()),
262                   x_ties_);
263
264    std::vector<double> vec(utility::pair_second_iterator(data.begin()),
265                            utility::pair_second_iterator(data.end()));
266    std::sort(vec.begin(), vec.end());
267    calculate_ties(vec.begin(), vec.end(), y_ties_);
268
269    /*
270                y1 < y2  y2 == y2  y2 > y2
271      x1 <  x2     C        eX        D
272      x1 == x2    eY       spare      -
273      x1 >  x2     -        -         -
274
275      We categorise pairs into five categories:
276      C: Concordant
277      D: Discordant
278      eX: extra X; Ys and only Ys are equal
279      eY: extra Y; Xs and only Xs are equal
280      spare: both Xs and Yy are equal
281
282      Due to symmetry reasons and because data container is sorted, we
283      can ignore lower part of the matrix above.
284     */
285
286    concordant_ = 0;
287    discordant_ = 0;
288    extraX_ = 0;
289    extraY_ = 0;
290
291    unsigned long int eY = 0;
292    // size of the current equal range, i.e., number of data points
293    // for X_i : X_j == X_i, Y_j == Y_i, j <= i including the current
294    // point
295    unsigned long int ties = 1; // because loop below skip first entry
296    utility::Ranking<double> Y;
297
298    // loop over data, which is sorted w.r.t. ::first
299    auto previous = data.cbegin();
300    assert(previous != data.cend());
301    Y.insert(previous->second);
302    auto it = std::next(previous);
303    while (it!=data.cend()) {
304      assert(previous->first <= it->first);
305      // X not equal
306      if (it->first != previous->first) {
307        eY = 0;
308        ties = 1;
309      }
310      // y also equal
311      else if (it->second == previous->second)
312        ++ties;
313      else { // x equal, y not equal
314        eY += ties;
315        ties = 1;
316      }
317
318      auto lower = Y.insert(it->second);
319      if (lower!=Y.begin() && *std::prev(lower) == it->second)
320        lower = Y.lower_bound(it->second);
321      // number of element in Y smaller than it->second
322      int n_smaller = Y.ranking(lower);
323      // number of element in Y equal to it->second
324      int n_equal = 1;
325      assert(lower != Y.cend());
326      auto upper = std::next(lower);
327      while (upper!=Y.cend() && *upper==*lower) {
328        ++upper;
329        ++n_equal;
330      }
331      size_t i = Y.size();
332
333      // n_smaller (y<yi) is the union of concordant (y<yi,x<xi)
334      // and eY (y<yi,x==xi)
335      int C = n_smaller - eY;
336
337      int eX =  n_equal - ties;
338
339      int D = i - (C + eX + eY + ties);
340
341      extraY_ += eY;
342      extraX_ += eX;
343      concordant_ += C;
344      discordant_ += D;
345      previous = it;
346      ++it;
347    }
348
349  }
350
351
352  long int Kendall::Pimpl::Count::count(void) const
353  {
354    return concordant_ - discordant_;
355  }
356
357
358  double Kendall::Pimpl::Count::score(void) const
359  {
360    double numerator = count();
361    double denominator = concordant_ + discordant_;
362    if (extraX_ || extraY_) {
363      denominator =
364        std::sqrt((denominator + extraX_)*(denominator + extraY_));
365    }
366    return numerator / denominator;
367  }
368
369
370  const Kendall::Pimpl::Count::Ties& Kendall::Pimpl::Count::x_ties(void) const
371  {
372    return x_ties_;
373  }
374
375
376  const Kendall::Pimpl::Count::Ties& Kendall::Pimpl::Count::y_ties(void) const
377  {
378    return y_ties_;
379  }
380
381
382  double Kendall::Pimpl::variance(void) const
383  {
384    /*
385      According to wikipedia,
386      z = k / sqrt(v)
387      is approximately standard normal
388      v = (v0 - vt - vu)/18 + v1 + v2
389      v0 = n(n-1)(2n+5)
390      vt = \sum t(t-1)(2t+5)
391      vu = \sum u(u-1)(2u+5)
392      v1 = \sum t(t-1)) * \sum u(u-1) / (2n(n-1))
393      v2 = sum t(t-1)(t-2) \sum u(u-1)(u-2) / (9n(n-1)(n-2))
394
395      where t is number of equal x values in group i and similarly u for
396      y.
397    */
398    double n = data_.size();
399    double v0 = n*(n-1)*(2*n+5);
400    double vt = 0;
401    double vu = 0;
402    double v1 = 0;
403    double v2 = 0;
404    assert(count_);
405    auto& x_ties = count_->x_ties();
406    auto& y_ties = count_->y_ties();
407    // all correction terms above are zero in absence of ties
408    bool x_have_ties = x_ties.have_ties();
409    bool y_have_ties = y_ties.have_ties();
410    if (x_have_ties || y_have_ties) {
411      if (x_have_ties)
412        vt = x_ties.v_correction();
413      if (y_have_ties) {
414        vu = y_ties.v_correction();
415        if (x_have_ties) {
416          v1 = x_ties.n_pairs() * (y_ties.n_pairs() / (2*n*(n-1)));
417          v2 = x_ties.n_triples();
418          if (v2)
419            v2 *= y_ties.n_triples() / (9*n*(n-1)*(n-2));
420        }
421      }
422    }
423    return (v0 - vt - vu)/18 + v1 + v2;
424  }
425
426
427  Kendall::Pimpl::Count::Ties::Ties(void)
428    : n_pairs_(0), n_triples_(0), v_correction_(0)
429  {}
430
431
432  void Kendall::Pimpl::Count::Ties::add(size_t n)
433  {
434    unsigned long int factor = n * (n-1);
435    n_pairs_ += factor;
436    n_triples_ += factor * (n-2);
437    v_correction_ += factor * (2*n+5);
438  }
439
440
441  Kendall::Pimpl::Pimpl(void)
442  {}
443
444
445  Kendall::Pimpl::Pimpl(const Pimpl& other)
446    : data_(other.data_)
447  {}
448
449
450  Kendall::Pimpl& Kendall::Pimpl::operator=(const Pimpl& rhs)
451  {
452    data_ = rhs.data_;
453    count_.reset();
454    return *this;
455  }
456
457
458  void Kendall::Pimpl::add(double x, double y)
459  {
460    data_.insert(std::make_pair(x, y));
461    count_.reset();
462  }
463
464
465  size_t Kendall::Pimpl::n(void) const
466  {
467    return data_.size();
468  }
469
470
471  double Kendall::Pimpl::p_approx(bool right) const
472  {
473    double k = count();
474    if (!right)
475      k = -k;
476    return gsl_cdf_gaussian_Q(k, std::sqrt(variance()));
477  }
478
479
480  double Kendall::Pimpl::p_exact(bool right, bool left) const
481  {
482    long int upper = 0;
483    long int lower = 0;
484    if (right) {
485      if (left) {
486        upper = std::max(count(), -count());
487        lower = -upper;
488      }
489      else {
490        upper = count();
491        lower = std::numeric_limits<long int>::min();
492      }
493    }
494    else {
495      assert(left && "left or right must be true");
496      upper = std::numeric_limits<long int>::max();
497      lower = count();
498    }
499
500    // create a copy of the data, sort it with respect to ::second and
501    // then iterate through the permutations of second while keeping
502    // first constant. It means we need to do one extra initial sort,
503    // but OTOH the permuted data is always almost sorted.
504    std::vector<std::pair<double,double>> data(data_.begin(), data_.end());
505    using utility::pair_second_iterator;
506    std::sort(pair_second_iterator(data.begin()),
507              pair_second_iterator(data.end()));
508    unsigned int n = 0;
509    unsigned int total = 0;
510    do {
511      std::multiset<std::pair<double,double>>
512                    dataset(data.begin(), data.end());
513      Count count(dataset);
514      if (count.count() <= lower || count.count() >= upper)
515        ++n;
516      ++total;
517    }
518    while (std::next_permutation(pair_second_iterator(data.begin()),
519                                 pair_second_iterator(data.end())));
520
521    return static_cast<double>(n)/static_cast<double>(total);
522  }
523
524
525  void Kendall::Pimpl::reset(void)
526  {
527    Pimpl tmp;
528    *this = tmp;
529  }
530
531
532  double Kendall::Pimpl::score(void) const
533  {
534    count();
535    assert(count_.get());
536    return count_->score();
537  }
538
539
540  long int Kendall::Pimpl::count(void) const
541  {
542    if (!count_)
543      // const_cast to allow lazy eval is more restrictive than
544      // making count_ mutable.
545      const_cast<Pimpl*>(this)->count_.reset(new Count(data_));
546    return count_->count();
547  }
548
549}}} // of namespace statistics, yat, and theplu
Note: See TracBrowser for help on using the repository browser.