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

Last change on this file since 303 was 303, checked in by Peter, 17 years ago

docs

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.3 KB
Line 
1// $Id: ROC.cc 303 2005-04-30 16:17:35Z peter $
2
3#include <c++_tools/statistics/ROC.h>
4
5#include <c++_tools/utility/stl_utility.h>
6#include <c++_tools/gslapi/vector.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 {
17namespace statistics { 
18
19  ROC::ROC(bool b) 
20    : Score(b), area_(-1), minimum_size_(10), nof_pos_(0), 
21      value_(std::vector<std::pair<double, double> >())
22     
23       
24  {
25  }
26
27  double ROC::get_p_approx(const double area) const
28  {
29    double x = area - 0.5;
30    // Not integrating from the middle of the bin, but from the inner edge.
31    if (x>0)
32      x -= 0.5/nof_pos_/(value_.size()-nof_pos_);
33    else if(x<0)
34      x += 0.5/nof_pos_/(value_.size()-nof_pos_);
35
36    double sigma = (std::sqrt( (value_.size()-nof_pos_)*nof_pos_*
37                               (value_.size()+1.0)/12 ) /
38                    ( value_.size() - nof_pos_ ) / nof_pos_ );
39    double p = gsl_cdf_gaussian_Q(x, sigma);
40   
41    return p;
42  }
43
44  double ROC::get_p_exact(const double block, const double nof_pos, 
45                          const double nof_neg) const
46  {
47    double p;
48    if (block <= 0.0)
49      p = 1.0;
50    else if (block > nof_neg*nof_pos)
51      p = 0.0;
52    else {
53      double p1 = get_p_exact(block-nof_neg, nof_pos-1, nof_neg);
54      double p2 = get_p_exact(block, nof_pos, nof_neg-1);
55      p = nof_pos/(nof_pos+nof_neg)*p1 + nof_neg/(nof_pos+nof_neg)*p2;
56    }
57    return p;
58  }
59
60  double ROC::p_value(void) const
61  {
62    if (weighted_)
63      return 1.0;
64    else if (nof_pos_ < minimum_size_ & value_.size()-nof_pos_ < minimum_size_)
65      return get_p_exact(area_*nof_pos_*(value_.size()-nof_pos_), 
66                          nof_pos_, value_.size()-nof_pos_);
67    else
68      return get_p_approx(area_);
69   
70  }
71
72  double ROC::score(const gslapi::vector& target, const gslapi::vector& data, 
73                    const std::vector<size_t>& train_set)
74  {
75    weighted_=false;
76    target_ = target;
77    data_ = data;
78    if (!train_set.size()){
79      train_set_.resize(0);
80      for (size_t i=0; i<target_.size(); i++)
81        train_set_.push_back(i); 
82    }
83    else
84      train_set_ = train_set;
85    sort();
86    area_ = 0;
87    for (size_t i=0; i<value_.size(); i++)
88      if (value_[i].first==1)
89        area_+=i;
90
91    // Normalizing the area to [0,1]
92    area_ = ( (area_-nof_pos_*(nof_pos_-1)/2 ) /
93              (nof_pos_*(value_.size()-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 gslapi::vector& target, const gslapi::vector& data, 
103                    const gslapi::vector& weight,
104                    const std::vector<size_t>& train_set)
105  {
106    weighted_=true;
107    target_ = target;
108    data_ = data;
109    weight_=weight;
110    if (!train_set.size()){
111      train_set_.resize(0);
112      for (size_t i=0; i<target_.size(); i++)
113        train_set_.push_back(i); 
114    }
115    else
116      train_set_ = train_set;
117    sort();
118    area_=0;
119    double max_area=0;
120    //Peter, use the sort to skip some ifs and loops
121    for (size_t i=0; i<value_.size(); i++){
122      if (target_(train_set_[i])==1){
123        for (size_t j=0; j<value_.size(); j++){
124          if (target_(train_set_[j])!=1){
125            if (data_(train_set_[i]) > data_(train_set_[j])){
126              area_+=weight_(train_set_[i])*weight_(train_set_[j]);
127            }
128            max_area+=weight_(train_set_[i])*weight_(train_set_[j]);
129          }
130        }
131      }
132    }
133    area_/=max_area;
134   
135    if (area_<0.5 && absolute_)
136      area_=1.0-area_;
137   
138    return area_;
139  }
140
141  void ROC::sort()
142  {
143    value_.resize(0);
144    nof_pos_=0;
145    for (unsigned int i=0; i<train_set_.size(); i++){
146      std::pair<double, double> tmp(target_(train_set_[i]), 
147                                    data_(train_set_[i]));
148      value_.push_back(tmp);
149      if (target_(train_set_[i])==1)
150        nof_pos_++;
151    }
152    std::sort(value_.begin(),value_.end(),
153         utility::pair_value_compare<double, double>());
154  }
155
156  gslapi::vector ROC::target(void) const
157  {
158    gslapi::vector target(value_.size());
159    for(size_t i=0; i<target.size(); ++i) 
160      target(i) = value_[i].first;
161    return target;
162  }
163
164  std::ostream& operator<<(std::ostream& s, const ROC& r)
165  {
166    s.setf( std::ios::dec );
167    s.precision(12);
168    double sens = 1;
169    double spec = 0;
170    gslapi::vector target = r.target();
171    size_t n = target.size();
172    double nof_pos = (target.sum()+n)/2;
173    for(size_t i=0; i<n-1; ++i) {
174      s << sens << "\t";
175      s << spec << "\n";
176      if (target(i)==1)
177        spec -= 1/(n-nof_pos);
178      else 
179        sens -= 1/nof_pos;
180    }
181    s << sens << "\t";
182    s << spec ;
183    return s;
184  }
185
186
187}} // of namespace statistics and namespace theplu
Note: See TracBrowser for help on using the repository browser.