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

Last change on this file since 4037 was 4037, checked in by Peter, 19 months ago

skeleton for new Ranking container. refs #710.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 9.8 KB
Line 
1// $Id: Kendall.cc 4037 2021-01-25 04:08:28Z 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      double variance(void) const;
84      class Ties
85      {
86      public:
87        Ties(void);
88        void add(size_t n);
89        bool have_ties(void) const { return n_pairs_; }
90        // \return \sum x * (x-1)
91        unsigned long int n_pairs(void) const { return n_pairs_;}
92        // \return \sum x * (x-1) * (x-2)
93        unsigned long int n_triples(void) const { return n_triples_; }
94        // \return \sum x * (x-1) * (2*x+5)
95        unsigned long int v_correction(void) const { return v_correction_; }
96      private:
97        unsigned long int n_pairs_;
98        unsigned long int n_triples_;
99        unsigned long int v_correction_;
100      };
101
102    private:
103      Ties x_ties_;
104      Ties y_ties_;
105
106      // # pairs such that (x_i-x_j)(y_i-y_j) > 0
107      long int concordant_;
108      // # pairs such that (x_i-x_j)(y_i-y_j) < 0
109      long int disconcordant_;
110      // # pairs such that x_i!=x_j && y_i==y_j
111      long int extraX_;
112      // # pairs such that x_i==x_j && y_i!=y_j
113      long int extraY_;
114      // # pairs such that x_i==x_j && y_i==y_j
115      long int spare_;
116
117      template<typename Iterator>
118      void calculate_ties(Iterator first, Iterator last, Ties& ties)
119      {
120        while (first != last) {
121          Iterator upper = first;
122          size_t n = 1;
123          ++upper;
124          while (upper!=last && *upper==*first) {
125            ++n;
126            ++upper;
127          }
128          ties.add(n);
129          first = upper;
130        }
131      }
132    };
133
134  public:
135    Pimpl(void);
136    Pimpl(const Pimpl& other);
137    Pimpl& operator=(const Pimpl& rhs);
138    void add(double x, double y);
139    size_t n(void) const;
140    /// \return one-sided p-value
141    double p_approx(bool right) const;
142    double p_exact(bool right, bool left) const;
143    void reset(void);
144    double score(void) const;
145  private:
146    // return # concordant pairs minus # discordant pairs
147    long int count(void) const;
148    // data always sort wrt first and then second (if first are equal)
149    std::multiset<std::pair<double, double> > data_;
150    // calculated in score(void)
151    boost::scoped_ptr<Count> count_;
152  };
153
154
155  // Kendall class
156  Kendall::Kendall(void)
157    : pimpl_(new Pimpl)
158  {
159  }
160
161
162  Kendall::Kendall(const Kendall& rhs)
163    : pimpl_(new Pimpl(*rhs.pimpl_))
164  {
165  }
166
167
168  Kendall::~Kendall(void)
169  {
170    delete pimpl_;
171  }
172
173
174  void Kendall::add(double x, double y)
175  {
176    pimpl_->add(x, y);
177  }
178
179
180  size_t Kendall::n(void) const
181  {
182    return pimpl_->n();
183  }
184
185
186  double Kendall::score(void) const
187  {
188    return pimpl_->score();
189  }
190
191
192  double Kendall::p_left(bool exact) const
193  {
194    if (!exact)
195      return pimpl_->p_approx(false);
196    return pimpl_->p_exact(false, true);
197  }
198
199
200  double Kendall::p_right(bool exact) const
201  {
202    if (!exact)
203      return pimpl_->p_approx(true);
204    return pimpl_->p_exact(true, false);
205  }
206
207
208  double Kendall::p_value(bool exact) const
209  {
210    if (exact)
211      return pimpl_->p_exact(true, true);
212    if (score()>0.0)
213      return 2*p_right(false);
214    return 2*p_left(false);
215  }
216
217
218  void Kendall::reset(void)
219  {
220    pimpl_->reset();
221  }
222
223
224  Kendall& Kendall::operator=(const Kendall& rhs)
225  {
226    if (&rhs == this)
227      return *this;
228
229    assert(pimpl_);
230    assert(rhs.pimpl_);
231    *pimpl_ = *rhs.pimpl_;
232    return *this;
233  }
234
235
236  Kendall::Pimpl::Count::Count(const std::multiset<std::pair<double,double>>& data)
237  {
238    // We follow 3 Algorithm SDTau for some-duplicate datasets in
239    // 'Fast Algorithms For The Calculation Of Kendall's Tau'
240    // by David Christen (Computational Statistics, March 2005)
241
242    // data is sorted w.r.t. ::first
243    calculate_ties(utility::pair_first_iterator(data.begin()),
244                   utility::pair_first_iterator(data.end()),
245                   x_ties_);
246    unsigned int d = 0;
247    unsigned int e = 0;
248    utility::Ranking<double> Y;
249
250    for (auto it=data.cbegin(); it!=data.end(); ++it) {
251      if (it != data.begin()) {
252        auto previous = std::prev(it);
253        if (it->first != previous->first) { // x not equal
254          d = 0;
255          e = 1;
256        }
257        else if (it->second == previous->second) // y also equal
258          ++e;
259        else { // x equal, y not equal
260          d += e;
261          e = 1;
262        }
263      }
264
265      auto lower = Y.lower_bound(it->second);
266      // number of element in Y smaller than it->second
267      int n_smaller = Y.ranking(lower);
268      // number of element in Y equal to it->second
269      int n_equal = 1;
270      auto upper = std::next(lower);
271      while (upper!=Y.cend() && *upper==*lower) {
272        ++upper;
273        ++n_equal;
274      }
275      Y.insert(lower, it->second);
276      size_t i = Y.size();
277
278      long int a = n_smaller - d;
279      long int b = n_equal - e;
280      long int c = i - (a+b+d+e-1);
281
282      extraY_ += d;
283      extraX_ += b;
284      concordant_ += a;
285      disconcordant_ += c;
286    }
287    assert(0);
288  }
289
290
291  long int Kendall::Pimpl::Count::count(void) const
292  {
293    return concordant_ - disconcordant_;
294  }
295
296
297  double Kendall::Pimpl::Count::score(void) const
298  {
299    double numerator = count();
300    double denominator = concordant_ + disconcordant_;
301    if (extraX_ || extraY_) {
302      denominator =
303        std::sqrt((denominator + extraX_)*(denominator + extraY_));
304    }
305    return numerator / denominator;
306  }
307
308
309  double Kendall::Pimpl::Count::variance(void) const
310  {
311    /*
312      According to wikipedia,
313      z = k / sqrt(v)
314      is approximately standard normal
315      v = (v0 - vt - vu)/18 + v1 + v2
316      v0 = n(n-1)(2n+5)
317      vt = \sum t(t-1)(2t+5)
318      vu = \sum u(u-1)(2u+5)
319      v1 = \sum t(t-1)) * \sum u(u-1) / (2n(n-1))
320      v2 = sum t(t-1)(t-2) \sum u(u-1)(u-2) / (9n(n-1)(n-2))
321
322      where t is number of equal values in group i and similarly u for
323      y.
324    */
325    double n = score();
326    double v0 = n*(n-1)*(2*n+5);
327    double vt = 0;
328    double vu = 0;
329    double v1 = 0;
330    double v2 = 0;
331    // all correction terms above are zero in absence of ties
332    bool x_have_ties = x_ties_.have_ties();
333    bool y_have_ties = y_ties_.have_ties();
334    if (x_have_ties || y_have_ties) {
335      if (x_have_ties)
336        vt = x_ties_.v_correction();
337      if (y_have_ties) {
338        vu = y_ties_.v_correction();
339        if (x_have_ties) {
340          v1 = x_ties_.n_pairs() * (y_ties_.n_pairs() / (2*n*(n-1)));
341          v2 = x_ties_.n_triples();
342          if (v2)
343            v2 *= y_ties_.n_triples() / (9*n*(n-1)*(n-2));
344        }
345      }
346    }
347    return (v0 - vt - vu)/18 + v1 + v2;
348  }
349
350
351  Kendall::Pimpl::Count::Ties::Ties(void)
352    : n_pairs_(0), n_triples_(0), v_correction_(0)
353  {}
354
355
356  void Kendall::Pimpl::Count::Ties::add(size_t n)
357  {
358    unsigned long int factor = n * (n-1);
359    n_pairs_ += factor;
360    n_triples_ += factor * (n-2);
361    v_correction_ += factor * (2*n+5);
362  }
363
364
365  Kendall::Pimpl::Pimpl(void)
366  {}
367
368
369  Kendall::Pimpl::Pimpl(const Pimpl& other)
370    : data_(other.data_)
371  {}
372
373
374  Kendall::Pimpl& Kendall::Pimpl::operator=(const Pimpl& rhs)
375  {
376    data_ = rhs.data_;
377    count_.reset();
378    return *this;
379  }
380
381
382  void Kendall::Pimpl::add(double x, double y)
383  {
384    data_.insert(std::make_pair(x, y));
385    count_.reset();
386  }
387
388
389  size_t Kendall::Pimpl::n(void) const
390  {
391    return data_.size();
392  }
393
394
395  double Kendall::Pimpl::p_approx(bool right) const
396  {
397    double k = count();
398    if (!right)
399      k = -k;
400    return gsl_cdf_gaussian_Q(k, std::sqrt(count_->variance()));
401  }
402
403
404  double Kendall::Pimpl::p_exact(bool right, bool left) const
405  {
406    assert(0);
407    return 0.0;
408    /*
409      long int upper = 0;
410      long int lower = 0;
411      if (right) {
412      if (left) {
413      upper = std::max(count(), -count());
414      lower = std::min(count(), -count());
415      }
416      else {
417      upper = count();
418      lower = std::numeric_limits<long int>::min();
419      }
420      }
421      else {
422      assert(left && "left or right must be true");
423      upper = std::numeric_limits<long int>::max();
424      lower = count();
425      }
426    */
427    /*
428      std::vector<double> x(x_);
429      std::sort(x.begin(), x.end());
430      unsigned int n = 0;
431      unsigned int total = 0;
432      do {
433      long int k = statistics::count(x.begin(), x.end(), y_.begin());
434      if (k>=upper || k<=lower)
435      ++n;
436      ++total;
437      }
438      while (std::next_permutation(x.begin(), x.end()));
439      return static_cast<double>(n)/static_cast<double>(total);
440    */
441    return 0;
442  }
443
444
445  void Kendall::Pimpl::reset(void)
446  {
447    Pimpl tmp;
448    *this = tmp;
449  }
450
451
452  double Kendall::Pimpl::score(void) const
453  {
454    count();
455    assert(count_.get());
456    return count_->score();
457  }
458
459
460  long int Kendall::Pimpl::count(void) const
461  {
462    if (!count_)
463      // const_cast to allow lazy eval is more restrictive than
464      // making count_ mutable.
465      const_cast<Pimpl*>(this)->count_.reset(new Count(data_));
466    return count_->count();
467  }
468
469}}} // of namespace statistics, yat, and theplu
Note: See TracBrowser for help on using the repository browser.