source: trunk/yat/statistics/ROC.cc @ 2595

Last change on this file since 2595 was 2595, checked in by Peter, 11 years ago

fixes #678. ROC p_value

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.9 KB
Line 
1// $Id: ROC.cc 2595 2011-10-30 03:27:20Z peter $
2
3/*
4  Copyright (C) 2004, 2005 Peter Johansson
5  Copyright (C) 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
6  Copyright (C) 2011 Peter Johansson
7
8  This file is part of the yat library, http://dev.thep.lu.se/yat
9
10  The yat library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU General Public License as
12  published by the Free Software Foundation; either version 3 of the
13  License, or (at your option) any later version.
14
15  The yat library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  General Public License for more details.
19
20  You should have received a copy of the GNU General Public License
21  along with yat. If not, see <http://www.gnu.org/licenses/>.
22*/
23
24#include "ROC.h"
25#include "AUC.h"
26
27#include <gsl/gsl_cdf.h>
28
29#include <algorithm>
30#include <cassert>
31#include <cmath>
32#include <utility>
33
34namespace theplu {
35namespace yat {
36namespace statistics { 
37
38  ROC::ROC(void) 
39    :minimum_size_(10)
40  {
41    reset();
42  }
43
44
45  void ROC::add(double x, bool target, double w)
46  {
47    if (!w)
48      return;
49    ROC::Map::value_type element(x, std::make_pair(target, w));
50    ROC::Map::iterator lower = multimap_.lower_bound(x);
51    if (lower!=multimap_.end() && lower->first == x)
52      has_ties_ = true;
53    multimap_.insert(lower, element);
54    if (target)
55      w_pos_+=w;
56    else
57      w_neg_+=w;
58    area_=-1;
59  }
60
61
62  double ROC::area(void)
63  {
64    if (area_==-1){
65      AUC auc(false);
66      area_=auc.score(multimap_);
67    }
68    return area_;
69  }
70
71
72  double ROC::get_p_approx(double x) const
73  {
74    // make x standard normal
75    x -= 0.5;
76    // Not integrating from the middle of the bin, but from the inner edge.
77    if (x>0)
78      x -= 0.5/(n_pos()*n_neg());
79    else if(x<0)
80      x += 0.5/(n_pos()*n_neg());
81    else
82      return 0.5;
83    double var = 1.0+n();
84    if (has_ties_) {
85      double correction = 0;
86      Map::const_iterator first = multimap_.begin();
87      Map::const_iterator last = multimap_.begin();
88      while (first!=multimap_.end()) {
89        size_t n = 0;
90        while (first->first == last->first) {
91          ++n;
92          ++last;
93        }
94        correction += n * (n-1) * (n+1);
95        first = last;
96      }
97      /*
98        mn(N+1)/12-[mn/(12N(N-1)) * sum(t(t-1)(t+1))] =
99        mn/12 [ N+1 - 1/(N(N-1)) * sum(t(t-1)(t+1)) ]
100      */
101      var -= correction/(n() * (n()-1));
102    }
103    var = var / (12*n_pos()*n_neg());
104    return gsl_cdf_gaussian_Q(x, std::sqrt(var));
105  }
106
107
108  unsigned int& ROC::minimum_size(void)
109  {
110    return minimum_size_;
111  }
112
113
114  const unsigned int& ROC::minimum_size(void) const
115  {
116    return minimum_size_;
117  }
118
119
120  double ROC::n(void) const
121  {
122    return n_pos()+n_neg();
123  }
124
125
126  double ROC::n_neg(void) const
127  {
128    return w_neg_;
129  }
130
131
132  double ROC::n_pos(void) const
133  {
134    return w_pos_;
135  }
136
137
138  double ROC::p_exact(double area) const
139  {
140    return p_exact_with_ties(multimap_.begin(), multimap_.end(),
141                             area*w_pos_*w_neg_, w_pos_, w_neg_);
142  }
143
144
145  double ROC::p_value() const
146  {
147    double area(area_);
148    if (area_==-1){
149      AUC auc(false);
150      area = auc.score(multimap_);
151    }
152    if (use_exact_method()) {
153      double p = 0;
154      double abs_area = std::max(area, 1-area);
155      p = p_exact(abs_area);
156      if (has_ties_) {
157        p += p_exact_with_ties(multimap_.rbegin(), multimap_.rend(),
158                               abs_area*w_pos_*w_neg_, w_pos_, w_neg_);
159      }
160      else
161        p *= 2.0;
162      // avoid double counting when area is 0.5
163      return std::min(p, 1.0);
164    }
165    return 2*get_p_approx(area);
166  }
167
168
169  double ROC::p_value_one_sided() const
170  {
171    double area(area_);
172    if (area_==-1){
173      AUC auc(false);
174      area = auc.score(multimap_);
175    }
176    if (use_exact_method())
177      return p_exact(area);
178    return get_p_approx(area);
179  }
180
181
182  void ROC::reset(void)
183  {
184    area_=-1;
185    has_ties_ = false;
186    w_pos_=0;
187    w_neg_=0;
188    multimap_.clear();
189  }
190
191
192  bool ROC::use_exact_method(void) const
193  {
194    return (n_pos() < minimum_size_) && (n_neg() < minimum_size_);
195  }
196
197}}} // of namespace statistics, yat, and theplu
Note: See TracBrowser for help on using the repository browser.