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

Last change on this file since 616 was 616, checked in by Jari Häkkinen, 15 years ago

Removed gslapi namespace and put the code into utility namespace.
Move #ifndef _header_ idiom to top of touched header files.
Removed unneccesary #includes, and added needed #includes.

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