source: trunk/src/ROC.cc @ 103

// $Id: ROC.cc 103 2004-06-15 14:29:15Z peter $

// System includes
#include <iostream>
//#include <algorithm>
//#include <utility>
//#include <vector>

// Thep C++ Tools
#include "ROC.h"
#include "stl_utility.h"
#include "vector.h"

namespace theplu {
namespace cpptools {  

  ROC::ROC(const gslapi::vector& target, const gslapi::vector& value) 
  
    : Score(), value_(), nof_pos_(0), minimum_size_(10), area_(-1)
    
  {
    sort(target, value);
  }

  ROC::ROC() 
    : Score(), value_(), nof_pos_(0), minimum_size_(10), area_(-1)
       
  {
  }

  double ROC::get_p_approx(const double area) const
  {
    double x = area - 0.5;
    // Not integrating from the middle of the bin, but from the inner edge.
    if (x>0)
      x -= 0.5/nof_pos_/(value_.size()-nof_pos_);
    else if(x<0)
      x += 0.5/nof_pos_/(value_.size()-nof_pos_);

    double sigma = (std::sqrt((value_.size()-nof_pos_)* nof_pos_ *
                        (value_.size()+1)/12) /
                    (value_.size() - nof_pos_ ) / nof_pos_);
    double p = gsl_cdf_gaussian_Q(x, sigma);
        
    return p;
  }

  double ROC::get_p_exact(const double block, const double nof_pos, 
                            const double nof_neg)
  {
    double p;
    if (block <= 0.0)
      p = 1.0;
    else if (block > nof_neg*nof_pos)
      p = 0.0;
    else {
      double p1 = get_p_exact(block-nof_neg, nof_pos-1, nof_neg);
      double p2 = get_p_exact(block, nof_pos, nof_neg-1);
      p = nof_pos/(nof_pos+nof_neg)*p1 + nof_neg/(nof_pos+nof_neg)*p2;
    }
    return p;
  }

  double ROC::p_value(void)
  {
    if (area_==-1)
      area_ = score();
    double p;
    if (nof_pos_ < minimum_size_ & value_.size()-nof_pos_ < minimum_size_)
      p = get_p_exact(area_*nof_pos_*(value_.size()-nof_pos_), 
                          nof_pos_, value_.size()-nof_pos_);
    else
    p = get_p_approx(area_);
    return p;
  }

  double ROC::score()
  {
    if (area_==-1){
      double area_=0;
      for (unsigned int i=0; i<value_.size(); i++)
        if (value_[i].first==1)
          area_+=i;
      // Normalizing the area to 0-1
      area_ = (area_/nof_pos_ - (nof_pos_ - 1)/2 )/(value_.size() - nof_pos_);
    }
    return area_;
  }

  double ROC::score(const gslapi::vector& target, const gslapi::vector& value)
  {
    sort(target, value);
    double area_=0;
    for (unsigned int i=0; i<value_.size(); i++)
      if (value_[i].first==1)
        area_+=i;
    // Normalizing the area to 0-1
    area_ = (area_/nof_pos_ - (nof_pos_ - 1)/2 )/(value_.size() - nof_pos_);

    return area_;
  }

  void ROC::sort(const gslapi::vector& target, const gslapi::vector& value)
  {
    for (unsigned int i=0; i<target.size(); i++){
      int targ=static_cast<int>(target(i));
      std::pair<int, double> tmp(targ, value(i));
      value_.push_back(tmp);
      if (targ==1)
        nof_pos_++;
    }
    std::sort(value_.begin(),value_.end(),
         pair_value_compare<int,double>());
  }

  gslapi::vector ROC::target(void) const
  {
    gslapi::vector target(value_.size());
    for(size_t i=0; i<target.size(); ++i) 
      target(i) = value_[i].first;
    return target;
  }

  std::ostream& operator<<(std::ostream& s, const ROC& r)
  {
    s.setf( std::ios::dec );
    s.precision(12);
    double sens = 1;
    double spec = 0;
    gslapi::vector target = r.target();
    size_t n = target.size();
    double nof_pos = (target.sum()+n)/2;
    for(size_t i=0; i<n-1; ++i) {
      s << sens << "\t";
      s << spec << std::endl;
      if (target(i)==1)
        spec -= 1/(n-nof_pos);
      else 
        sens -= 1/nof_pos;
    }
    s << sens << "\t";
    s << spec << std::endl;
    return s;
  }


}} // of namespace cpptools and namespace theplu