1 | // $Id: ROC.cc 469 2005-12-19 14:58:29Z peter $ |
---|
2 | |
---|
3 | #include <c++_tools/statistics/ROC.h> |
---|
4 | |
---|
5 | #include <c++_tools/utility/stl_utility.h> |
---|
6 | #include <c++_tools/classifier/VectorAbstract.h> |
---|
7 | |
---|
8 | #include <gsl/gsl_cdf.h> |
---|
9 | |
---|
10 | #include <cmath> |
---|
11 | #include <iostream> |
---|
12 | #include <utility> |
---|
13 | #include <vector> |
---|
14 | |
---|
15 | |
---|
16 | namespace theplu { |
---|
17 | class classifier::VectorAbstract; |
---|
18 | namespace statistics { |
---|
19 | |
---|
20 | ROC::ROC(bool b) |
---|
21 | : Score(b), area_(0.5), minimum_size_(10), nof_pos_(0) |
---|
22 | { |
---|
23 | } |
---|
24 | |
---|
25 | double ROC::get_p_approx(const double area) const |
---|
26 | { |
---|
27 | double x = area - 0.5; |
---|
28 | // Not integrating from the middle of the bin, but from the inner edge. |
---|
29 | if (x>0) |
---|
30 | x -= 0.5/nof_pos_/(n()-nof_pos_); |
---|
31 | else if(x<0) |
---|
32 | x += 0.5/nof_pos_/(n()-nof_pos_); |
---|
33 | |
---|
34 | double sigma = (std::sqrt( (n()-nof_pos_)*nof_pos_* |
---|
35 | (n()+1.0)/12 ) / |
---|
36 | ( n() - nof_pos_ ) / nof_pos_ ); |
---|
37 | double p = gsl_cdf_gaussian_Q(x, sigma); |
---|
38 | |
---|
39 | return p; |
---|
40 | } |
---|
41 | |
---|
42 | double ROC::get_p_exact(const double block, const double nof_pos, |
---|
43 | const double nof_neg) const |
---|
44 | { |
---|
45 | double p; |
---|
46 | if (block <= 0.0) |
---|
47 | p = 1.0; |
---|
48 | else if (block > nof_neg*nof_pos) |
---|
49 | p = 0.0; |
---|
50 | else { |
---|
51 | double p1 = get_p_exact(block-nof_neg, nof_pos-1, nof_neg); |
---|
52 | double p2 = get_p_exact(block, nof_pos, nof_neg-1); |
---|
53 | p = nof_pos/(nof_pos+nof_neg)*p1 + nof_neg/(nof_pos+nof_neg)*p2; |
---|
54 | } |
---|
55 | return p; |
---|
56 | } |
---|
57 | |
---|
58 | double ROC::p_value(void) const |
---|
59 | { |
---|
60 | if (weighted_) |
---|
61 | return 1.0; |
---|
62 | else if (nof_pos_ < minimum_size_ & n()-nof_pos_ < minimum_size_) |
---|
63 | return get_p_exact(area_*nof_pos_*(n()-nof_pos_), |
---|
64 | nof_pos_, n()-nof_pos_); |
---|
65 | else |
---|
66 | return get_p_approx(area_); |
---|
67 | |
---|
68 | } |
---|
69 | |
---|
70 | double ROC::score(const classifier::Target& target, |
---|
71 | const classifier::VectorAbstract& value) |
---|
72 | { |
---|
73 | weighted_=false; |
---|
74 | |
---|
75 | vec_pair_.clear(); |
---|
76 | vec_pair_.reserve(target.size()); |
---|
77 | for (unsigned int i=0; i<target.size(); i++) |
---|
78 | vec_pair_.push_back(std::make_pair(target(i),value(i))); |
---|
79 | std::sort(vec_pair_.begin(),vec_pair_.end(), |
---|
80 | utility::pair_value_compare<int, double>()); |
---|
81 | |
---|
82 | area_ = 0; |
---|
83 | nof_pos_=0; |
---|
84 | for (size_t i=0; i<n(); i++){ |
---|
85 | if (class_one(target(i))){ |
---|
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 | double ROC::score(const classifier::Target& target, |
---|
103 | const classifier::VectorAbstract& value, |
---|
104 | const classifier::VectorAbstract& 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 | vec_pair_.push_back(std::make_pair(target(i),value(i))); |
---|
112 | std::sort(vec_pair_.begin(),vec_pair_.end(), |
---|
113 | utility::pair_value_compare<int, double>()); |
---|
114 | |
---|
115 | area_=0; |
---|
116 | nof_pos_=0; |
---|
117 | double max_area=0; |
---|
118 | |
---|
119 | for (size_t i=0; i<n(); i++){ |
---|
120 | if (class_one(target(i))){ |
---|
121 | for (size_t j=0; j<n(); j++){ |
---|
122 | if (!class_one(target(j))){ |
---|
123 | if (value(i) > value(j)){ |
---|
124 | area_+=weight(i)*weight(j); |
---|
125 | } |
---|
126 | max_area+=weight(i)*weight(j); |
---|
127 | } |
---|
128 | } |
---|
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 | int ROC::target(const size_t i) const |
---|
140 | { |
---|
141 | return vec_pair_[i].first; |
---|
142 | } |
---|
143 | |
---|
144 | std::ostream& operator<<(std::ostream& s, const ROC& r) |
---|
145 | { |
---|
146 | s.setf( std::ios::dec ); |
---|
147 | s.precision(12); |
---|
148 | double sens = 1; |
---|
149 | double spec = 0; |
---|
150 | size_t n = r.n(); |
---|
151 | double nof_pos = r.n_pos(); |
---|
152 | for(size_t i=0; i<n-1; ++i) { |
---|
153 | s << sens << "\t"; |
---|
154 | s << spec << "\n"; |
---|
155 | if (r.class_one(r.target(i))) |
---|
156 | spec -= 1/(n-nof_pos); |
---|
157 | else |
---|
158 | sens -= 1/nof_pos; |
---|
159 | } |
---|
160 | s << sens << "\t"; |
---|
161 | s << spec ; |
---|
162 | return s; |
---|
163 | } |
---|
164 | |
---|
165 | |
---|
166 | }} // of namespace statistics and namespace theplu |
---|