source: trunk/src/ROC.cc @ 112

Last change on this file since 112 was 112, checked in by Peter, 18 years ago

added the choice to not use all data points but just the train_set

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.7 KB
Line 
1// $Id: ROC.cc 112 2004-07-07 10:23:44Z peter $
2
3// System includes
4#include <iostream>
5//#include <algorithm>
6#include <utility>
7#include <vector>
8#include <cmath>
9
10// Thep C++ Tools
11#include "ROC.h"
12#include "stl_utility.h"
13#include "vector.h"
14
15namespace theplu {
16namespace cpptools { 
17
18  ROC::ROC(const gslapi::vector& target, const gslapi::vector& data, 
19           const std::vector<size_t>& train_set) 
20 
21    : Score(), area_(-1), data_(data), minimum_size_(10), nof_pos_(0), 
22      target_(target), train_set_(train_set), 
23      value_(std::vector<std::pair<double, double> >())
24   
25  {
26    if (!train_set_.size())
27      for (size_t i=0; i<target_.size(); i++)
28        train_set_.push_back(i); 
29    sort();
30  }
31
32  ROC::ROC() 
33    : Score(), area_(-1), data_(), minimum_size_(10), nof_pos_(0), target_(), 
34      train_set_(std::vector<size_t>()), 
35      value_(std::vector<std::pair<double, double> >())
36       
37  {
38  }
39
40  double ROC::get_p_approx(const double area) const
41  {
42    double x = area - 0.5;
43    // Not integrating from the middle of the bin, but from the inner edge.
44    if (x>0)
45      x -= 0.5/nof_pos_/(value_.size()-nof_pos_);
46    else if(x<0)
47      x += 0.5/nof_pos_/(value_.size()-nof_pos_);
48
49    double sigma = (std::sqrt( (value_.size()-nof_pos_)*nof_pos_*
50                               (value_.size()+1.0)/12 ) /
51                    ( value_.size() - nof_pos_ ) / nof_pos_ );
52    double p = gsl_cdf_gaussian_Q(x, sigma);
53       
54    return p;
55  }
56
57  double ROC::get_p_exact(const double block, const double nof_pos, 
58                            const double nof_neg)
59  {
60    double p;
61    if (block <= 0.0)
62      p = 1.0;
63    else if (block > nof_neg*nof_pos)
64      p = 0.0;
65    else {
66      double p1 = get_p_exact(block-nof_neg, nof_pos-1, nof_neg);
67      double p2 = get_p_exact(block, nof_pos, nof_neg-1);
68      p = nof_pos/(nof_pos+nof_neg)*p1 + nof_neg/(nof_pos+nof_neg)*p2;
69    }
70    return p;
71  }
72
73  double ROC::p_value(void)
74  {
75    if (area_==-1)
76      area_ = score();
77    double p;
78    if (nof_pos_ < minimum_size_ & value_.size()-nof_pos_ < minimum_size_)
79      p = get_p_exact(area_*nof_pos_*(value_.size()-nof_pos_), 
80                          nof_pos_, value_.size()-nof_pos_);
81    else
82    p = get_p_approx(area_);
83    return p;
84  }
85
86  double ROC::score()
87  {
88    if (area_==-1){
89      double area_ = 0;
90      for (unsigned int i=0; i<value_.size(); i++)
91        if (value_[i].first==1)
92          area_+=i;
93      // Normalizing the area to 0-1
94      area_ = (area_/nof_pos_ - (nof_pos_ - 1)/2 )/(value_.size() - nof_pos_);
95    }
96    return area_;
97  }
98
99  double ROC::score(const gslapi::vector& target, const gslapi::vector& data, 
100                    const std::vector<size_t>& train_set)
101  {
102    target_ = target;
103    data_ = data;
104    if (!train_set.size()){
105      train_set_.resize(0);
106      for (size_t i=0; i<target_.size(); i++)
107        train_set_.push_back(i); 
108    }
109    else
110      train_set_ = train_set;
111    sort();
112    area_ = score();   
113    return area_;
114  }
115
116  void ROC::sort()
117  {
118    value_.resize(0);
119    for (unsigned int i=0; i<train_set_.size(); i++){
120      std::pair<double, double> tmp(target_(train_set_[i]), data_(train_set_[i]));
121      value_.push_back(tmp);
122      if (target_(train_set_[i])==1)
123        nof_pos_++;
124    }
125    std::sort(value_.begin(),value_.end(),
126         pair_value_compare<double, double>());
127  }
128
129  gslapi::vector ROC::target(void) const
130  {
131    gslapi::vector target(value_.size());
132    for(size_t i=0; i<target.size(); ++i) 
133      target(i) = value_[i].first;
134    return target;
135  }
136
137  std::ostream& operator<<(std::ostream& s, const ROC& r)
138  {
139    s.setf( std::ios::dec );
140    s.precision(12);
141    double sens = 1;
142    double spec = 0;
143    gslapi::vector target = r.target();
144    size_t n = target.size();
145    double nof_pos = (target.sum()+n)/2;
146    for(size_t i=0; i<n-1; ++i) {
147      s << sens << "\t";
148      s << spec << "\n";
149      if (target(i)==1)
150        spec -= 1/(n-nof_pos);
151      else 
152        sens -= 1/nof_pos;
153    }
154    s << sens << "\t";
155    s << spec ;
156    return s;
157  }
158
159
160}} // of namespace cpptools and namespace theplu
Note: See TracBrowser for help on using the repository browser.