source: branches/peters_vector/lib/statistics/ROC.cc @ 470

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

compiling 3 tests fail

  • 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 470 2005-12-19 19:46:50Z 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 <iostream>
11#include <utility>
12#include <vector>
13
14
15namespace theplu {
16  class gslapi::vector;
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 gslapi::vector& value)
71  {
72    weighted_=false;
73
74    vec_pair_.clear();
75    vec_pair_.reserve(target.size());
76    for (unsigned int i=0; i<target.size(); i++)
77      vec_pair_.push_back(std::make_pair(target(i),value(i)));
78    std::sort(vec_pair_.begin(),vec_pair_.end(),
79              utility::pair_value_compare<int, double>());
80
81    area_ = 0;
82    nof_pos_=0;
83    for (size_t i=0; i<n(); i++){
84      if (class_one(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  double ROC::score(const classifier::Target& target, 
102                    const gslapi::vector& value, 
103                    const gslapi::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      vec_pair_.push_back(std::make_pair(target(i),value(i)));
111    std::sort(vec_pair_.begin(),vec_pair_.end(),
112              utility::pair_value_compare<int, double>());
113
114    area_=0;
115    nof_pos_=0;
116    double max_area=0;
117
118    for (size_t i=0; i<n(); i++){
119      if (class_one(vec_pair_[i].first)){
120        for (size_t j=0; j<n(); j++){
121          if (!class_one(vec_pair_[j].first)){
122            if (value(i) > value(j)){
123              area_+=weight(i)*weight(j);
124            }
125            max_area+=weight(i)*weight(j);
126          }
127        }
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  int 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.class_one(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.