source: trunk/lib/statistics/ROC.cc @ 509

Last change on this file since 509 was 509, checked in by Peter, 16 years ago

added test for target
redesign crossSplitter
added two class function in Target

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