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

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

non compiling checking before revision after design meeting

  • 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 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
16namespace theplu {
17  class classifier::VectorAbstract;
18namespace 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
Note: See TracBrowser for help on using the repository browser.