source: trunk/src/ROC.cc @ 110

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

fixed output operator

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