source: trunk/c++_tools/statistics/ROC.cc @ 623

Last change on this file since 623 was 623, checked in by Peter, 15 years ago

fixes #112 and refs #123 added overloaded function score taking Target and DataLookupWeighted1D, which is needed for InputRanker?.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.5 KB
Line 
1// $Id: ROC.cc 623 2006-09-05 02:13:12Z peter $
2
3#include "c++_tools/statistics/ROC.h"
4
5#include "c++_tools/classifier/DataLookupWeighted1D.h"
6#include "c++_tools/utility/stl_utility.h"
7#include "c++_tools/utility/vector.h"
8
9#include <gsl/gsl_cdf.h>
10
11#include <cmath>
12#include <utility>
13#include <vector>
14
15
16namespace theplu {
17namespace statistics { 
18
19  ROC::ROC(bool b) 
20    : Score(b), area_(0.5), minimum_size_(10), nof_pos_(0) 
21  {
22  }
23
24  double ROC::get_p_approx(const double area) const
25  {
26    double x = area - 0.5;
27    // Not integrating from the middle of the bin, but from the inner edge.
28    if (x>0)
29      x -= 0.5/nof_pos_/(n()-nof_pos_);
30    else if(x<0)
31      x += 0.5/nof_pos_/(n()-nof_pos_);
32
33    double sigma = (std::sqrt( (n()-nof_pos_)*nof_pos_*
34                               (n()+1.0)/12 ) /
35                    ( n() - nof_pos_ ) / nof_pos_ );
36    double p = gsl_cdf_gaussian_Q(x, sigma);
37   
38    return p;
39  }
40
41  double ROC::get_p_exact(const double block, const double nof_pos, 
42                          const double nof_neg) const
43  {
44    double p;
45    if (block <= 0.0)
46      p = 1.0;
47    else if (block > nof_neg*nof_pos)
48      p = 0.0;
49    else {
50      double p1 = get_p_exact(block-nof_neg, nof_pos-1, nof_neg);
51      double p2 = get_p_exact(block, nof_pos, nof_neg-1);
52      p = nof_pos/(nof_pos+nof_neg)*p1 + nof_neg/(nof_pos+nof_neg)*p2;
53    }
54    return p;
55  }
56
57  double ROC::p_value(void) const
58  {
59    if (weighted_)
60      return 1.0;
61    else if (nof_pos_ < minimum_size_ & n()-nof_pos_ < minimum_size_)
62      return get_p_exact(area_*nof_pos_*(n()-nof_pos_), 
63                          nof_pos_, n()-nof_pos_);
64    else
65      return get_p_approx(area_);
66   
67  }
68
69  double ROC::score(const classifier::Target& target, 
70                    const utility::vector& value)
71  {
72    assert(target.size()==value.size());
73    weighted_=false;
74
75    vec_pair_.clear();
76    vec_pair_.reserve(target.size());
77    for (size_t i=0; i<target.size(); i++)
78      vec_pair_.push_back(std::make_pair(target.binary(i),value(i)));
79
80    std::sort(vec_pair_.begin(),vec_pair_.end(),
81              utility::pair_value_compare<bool, double>());
82    area_ = 0;
83    nof_pos_=0;
84    for (size_t i=0; i<n(); i++){
85      if (vec_pair_[i].first){
86        area_+=i;
87        nof_pos_++;
88      }
89    }
90
91    // Normalizing the area to [0,1]
92    area_ = ( (area_-nof_pos_*(nof_pos_-1)/2 ) /
93              (nof_pos_*(n()-nof_pos_)) );
94   
95    //Returning score larger 0.5 that you get by random
96    if (area_<0.5 && absolute_)
97      area_=1.0-area_;
98
99    return area_;
100  }
101
102
103  // Peter, should be possible to do this in NlogN
104  double ROC::score(const classifier::Target& target, 
105                    const classifier::DataLookupWeighted1D& value)
106  {
107    weighted_=true;
108
109    vec_pair_.clear();
110    vec_pair_.reserve(target.size());
111    for (unsigned int i=0; i<target.size(); i++)
112      if (value.weight(i))
113        vec_pair_.push_back(std::make_pair(target.binary(i),value.data(i)));
114
115    std::sort(vec_pair_.begin(),vec_pair_.end(),
116              utility::pair_value_compare<int, double>());
117
118    area_=0;
119    nof_pos_=0;
120    double max_area=0;
121
122    for (size_t i=0; i<n(); i++)
123      if (target.binary(i))
124        for (size_t j=0; j<n(); j++)
125          if (!target.binary(j)){
126            if (value.data(i)>value.data(j))
127              area_+=value.weight(i)*value.weight(j);
128            max_area+=value.weight(i)*value.weight(j);
129          }
130   
131    area_/=max_area;
132   
133    if (area_<0.5 && absolute_)
134      area_=1.0-area_;
135   
136    return area_;
137  }
138
139
140  // Peter, should be possible to do this in NlogN
141  double ROC::score(const classifier::Target& target, 
142                    const utility::vector& value, 
143                    const utility::vector& weight)
144  {
145    weighted_=true;
146
147    vec_pair_.clear();
148    vec_pair_.reserve(target.size());
149    for (unsigned int i=0; i<target.size(); i++)
150      if (weight(i))
151        vec_pair_.push_back(std::make_pair(target.binary(i),value(i)));
152
153    std::sort(vec_pair_.begin(),vec_pair_.end(),
154              utility::pair_value_compare<int, double>());
155
156    area_=0;
157    nof_pos_=0;
158    double max_area=0;
159
160    for (size_t i=0; i<n(); i++)
161      if (target.binary(i))
162        for (size_t j=0; j<n(); j++)
163          if (!target.binary(j)){
164            if (value(i)>value(j))
165              area_+=weight(i)*weight(j);
166            max_area+=weight(i)*weight(j);
167          }
168   
169    area_/=max_area;
170   
171    if (area_<0.5 && absolute_)
172      area_=1.0-area_;
173   
174    return area_;
175  }
176
177  bool ROC::target(const size_t i) const
178  {
179    return vec_pair_[i].first;
180  }
181
182  std::ostream& operator<<(std::ostream& s, const ROC& r)
183  {
184    s.setf( std::ios::dec );
185    s.precision(12);
186    double sens = 1;
187    double spec = 0;
188    size_t n = r.n();
189    double nof_pos = r.n_pos();
190    for(size_t i=0; i<n-1; ++i) {
191      s << sens << "\t";
192      s << spec << "\n";
193      if (r.target(i))
194        spec -= 1/(n-nof_pos);
195      else 
196        sens -= 1/nof_pos;
197    }
198    s << sens << "\t";
199    s << spec ;
200    return s;
201  }
202
203
204}} // of namespace statistics and namespace theplu
Note: See TracBrowser for help on using the repository browser.