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

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

made Score classes ignoring data points with weight zero

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.7 KB
Line 
1// $Id: ROC.cc 488 2006-01-04 17:30:45Z 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    weighted_=false;
72
73    vec_pair_.clear();
74    vec_pair_.reserve(target.size());
75    for (unsigned int i=0; i<target.size(); i++){
76      int target_tmp = target(i);
77      double value_tmp = value(i);
78      vec_pair_.push_back(std::make_pair(target_tmp,value_tmp));
79      //      vec_pair_.push_back(std::make_pair(target(i),value(i)));
80    }
81    std::sort(vec_pair_.begin(),vec_pair_.end(),
82              utility::pair_value_compare<int, double>());
83
84
85    area_ = 0;
86    nof_pos_=0;
87    for (size_t i=0; i<n(); i++){
88      if (class_one(vec_pair_[i].first)){
89        area_+=i;
90        nof_pos_++;
91      }
92    }
93
94    // Normalizing the area to [0,1]
95    area_ = ( (area_-nof_pos_*(nof_pos_-1)/2 ) /
96              (nof_pos_*(n()-nof_pos_)) );
97   
98    //Returning score larger 0.5 that you get by random
99    if (area_<0.5 && absolute_)
100      area_=1.0-area_;
101
102    return area_;
103  }
104
105  // Peter, should be possible to do this in NlogN
106  double ROC::score(const classifier::Target& target, 
107                    const gslapi::vector& value, 
108                    const gslapi::vector& weight)
109  {
110    weighted_=true;
111
112    vec_pair_.clear();
113    vec_pair_.reserve(target.size());
114    for (unsigned int i=0; i<target.size(); i++)
115      if (weight(i))
116        vec_pair_.push_back(std::make_pair(target(i),value(i)));
117
118    std::sort(vec_pair_.begin(),vec_pair_.end(),
119              utility::pair_value_compare<int, double>());
120
121    area_=0;
122    nof_pos_=0;
123    double max_area=0;
124
125    for (size_t i=0; i<n(); i++)
126      if (class_one(target(i)))
127        for (size_t j=0; j<n(); j++)
128          if (!class_one(target(j))){
129            if (value(i)>value(j))
130              area_+=weight(i)*weight(j);
131            max_area+=weight(i)*weight(j);
132          }
133   
134    area_/=max_area;
135   
136    if (area_<0.5 && absolute_)
137      area_=1.0-area_;
138   
139    return area_;
140  }
141
142  int ROC::target(const size_t i) const
143  {
144    return vec_pair_[i].first;
145  }
146
147  std::ostream& operator<<(std::ostream& s, const ROC& r)
148  {
149    s.setf( std::ios::dec );
150    s.precision(12);
151    double sens = 1;
152    double spec = 0;
153    size_t n = r.n();
154    double nof_pos = r.n_pos();
155    for(size_t i=0; i<n-1; ++i) {
156      s << sens << "\t";
157      s << spec << "\n";
158      if (r.class_one(r.target(i)))
159        spec -= 1/(n-nof_pos);
160      else 
161        sens -= 1/nof_pos;
162    }
163    s << sens << "\t";
164    s << spec ;
165    return s;
166  }
167
168
169}} // of namespace statistics and namespace theplu
Note: See TracBrowser for help on using the repository browser.