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

// $Id: ROC.cc 303 2005-04-30 16:17:35Z peter $

#include <c++_tools/statistics/ROC.h>

#include <c++_tools/utility/stl_utility.h>
#include <c++_tools/gslapi/vector.h>

#include <gsl/gsl_cdf.h>

#include <cmath>
#include <iostream>
#include <utility>
#include <vector>


namespace theplu {
namespace statistics {  

  ROC::ROC(bool b) 
    : Score(b), area_(-1), minimum_size_(10), nof_pos_(0),  
      value_(std::vector<std::pair<double, double> >())
      
       
  {
  }

  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.0)/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) const
  {
    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) const
  {
    if (weighted_)
      return 1.0;
    else if (nof_pos_ < minimum_size_ & value_.size()-nof_pos_ < minimum_size_)
      return get_p_exact(area_*nof_pos_*(value_.size()-nof_pos_), 
                          nof_pos_, value_.size()-nof_pos_);
    else
      return get_p_approx(area_);
   
  }

  double ROC::score(const gslapi::vector& target, const gslapi::vector& data, 
                    const std::vector<size_t>& train_set)
  {
    weighted_=false;
    target_ = target;
    data_ = data;
    if (!train_set.size()){
      train_set_.resize(0);
      for (size_t i=0; i<target_.size(); i++)
        train_set_.push_back(i);  
    }
    else
      train_set_ = train_set;
    sort();
    area_ = 0;
    for (size_t 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 ) /
              (nof_pos_*(value_.size()-nof_pos_)) );
    
    //Returning score larger 0.5 that you get by random
    if (area_<0.5 && absolute_)
      area_=1.0-area_;
    
    return area_;
  }

  double ROC::score(const gslapi::vector& target, const gslapi::vector& data, 
                    const gslapi::vector& weight,
                    const std::vector<size_t>& train_set)
  {
    weighted_=true;
    target_ = target;
    data_ = data;
    weight_=weight;
    if (!train_set.size()){
      train_set_.resize(0);
      for (size_t i=0; i<target_.size(); i++)
        train_set_.push_back(i);  
    }
    else
      train_set_ = train_set;
    sort();
    area_=0;
    double max_area=0;
    //Peter, use the sort to skip some ifs and loops
    for (size_t i=0; i<value_.size(); i++){
      if (target_(train_set_[i])==1){
        for (size_t j=0; j<value_.size(); j++){
          if (target_(train_set_[j])!=1){
            if (data_(train_set_[i]) > data_(train_set_[j])){
              area_+=weight_(train_set_[i])*weight_(train_set_[j]);
            }
            max_area+=weight_(train_set_[i])*weight_(train_set_[j]);
          }
        }
      }
    }
    area_/=max_area;
    
    if (area_<0.5 && absolute_)
      area_=1.0-area_;
    
    return area_;
  }

  void ROC::sort()
  {
    value_.resize(0);
    nof_pos_=0;
    for (unsigned int i=0; i<train_set_.size(); i++){
      std::pair<double, double> tmp(target_(train_set_[i]), 
                                    data_(train_set_[i]));
      value_.push_back(tmp);
      if (target_(train_set_[i])==1)
        nof_pos_++;
    }
    std::sort(value_.begin(),value_.end(),
         utility::pair_value_compare<double, 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 << "\n";
      if (target(i)==1)
        spec -= 1/(n-nof_pos);
      else 
        sens -= 1/nof_pos;
    }
    s << sens << "\t";
    s << spec ;
    return s;
  }


}} // of namespace statistics and namespace theplu