source: trunk/yat/statistics/ROC.cc @ 781

// $Id: ROC.cc 781 2007-03-05 19:44:03Z peter $

/*
  Copyright (C) The authors contributing to this file.

  This file is part of the yat library, http://lev.thep.lu.se/trac/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "ROC.h"
#include "yat/classifier/DataLookupWeighted1D.h"
#include "yat/classifier/Target.h"
#include "yat/utility/stl_utility.h"
#include "yat/utility/vector.h"

#include <gsl/gsl_cdf.h>

#include <cassert>
#include <cmath>
#include <utility>
#include <vector>

namespace theplu {
namespace yat {
namespace statistics {  

  ROC::ROC(void) 
    : area_(0.5), minimum_size_(10), nof_pos_(0)  
  {
  }

  ROC::~ROC(void)
  {
  }

  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_/(n()-nof_pos_);
    else if(x<0)
      x += 0.5/nof_pos_/(n()-nof_pos_);

    double sigma = (std::sqrt( (n()-nof_pos_)*nof_pos_*
                               (n()+1.0)/12 ) /
                    ( n() - 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;
  }

  u_int& ROC::minimum_size(void)
  {
    return minimum_size_;
  }

  size_t ROC::n(void) const
  {
    return vec_pair_.size();
  }

  size_t ROC::n_pos(void) const
  {
    return nof_pos_;
  }

  double ROC::p_value(void) const
  {
    if (weighted_)
      return 1.0;
    else if (nof_pos_ < minimum_size_ & n()-nof_pos_ < minimum_size_)
      return get_p_exact(area_*nof_pos_*(n()-nof_pos_), 
                          nof_pos_, n()-nof_pos_);
    else
      return get_p_approx(area_);
   
  }

  double ROC::score(const classifier::Target& target, 
                    const utility::vector& value)
  {
    assert(target.size()==value.size());
    weighted_=false;

    vec_pair_.clear();
    vec_pair_.reserve(target.size());
    for (size_t i=0; i<target.size(); i++)
      vec_pair_.push_back(std::make_pair(target.binary(i),value(i)));

    std::sort(vec_pair_.begin(),vec_pair_.end(),
              utility::pair_value_compare<bool, double>());
    area_ = 0;
    nof_pos_=0;
    for (size_t i=0; i<n(); i++){
      if (vec_pair_[i].first){
        area_+=i;
        nof_pos_++;
      }
    }

    // Normalizing the area to [0,1]
    area_ = ( (area_-nof_pos_*(nof_pos_-1)/2 ) /
              (nof_pos_*(n()-nof_pos_)) );
    
    return area_;
  }


  // Peter, should be possible to do this in NlogN
  double ROC::score(const classifier::Target& target, 
                    const classifier::DataLookupWeighted1D& value)
  {
    weighted_=true;

    vec_pair_.clear();
    vec_pair_.reserve(target.size());
    for (unsigned int i=0; i<target.size(); i++)
      if (value.weight(i))
        vec_pair_.push_back(std::make_pair(target.binary(i),value.data(i)));

    std::sort(vec_pair_.begin(),vec_pair_.end(),
              utility::pair_value_compare<int, double>());

    area_=0;
    nof_pos_=0;
    double max_area=0;

    for (size_t i=0; i<n(); i++)
      if (target.binary(i))
        for (size_t j=0; j<n(); j++)
          if (!target.binary(j)){
            if (value.data(i)>value.data(j))
              area_+=value.weight(i)*value.weight(j);
            max_area+=value.weight(i)*value.weight(j);
          }
    
    area_/=max_area;
    
    return area_;
  }


  // Peter, should be possible to do this in NlogN
  double ROC::score(const classifier::Target& target, 
                    const utility::vector& value, 
                    const utility::vector& weight)
  {
    weighted_=true;

    vec_pair_.clear();
    vec_pair_.reserve(target.size());
    for (unsigned int i=0; i<target.size(); i++)
      if (weight(i))
        vec_pair_.push_back(std::make_pair(target.binary(i),value(i)));

    std::sort(vec_pair_.begin(),vec_pair_.end(),
              utility::pair_value_compare<int, double>());

    area_=0;
    nof_pos_=0;
    double max_area=0;

    for (size_t i=0; i<n(); i++)
      if (target.binary(i))
        for (size_t j=0; j<n(); j++)
          if (!target.binary(j)){
            if (value(i)>value(j))
              area_+=weight(i)*weight(j);
            max_area+=weight(i)*weight(j);
          }
    
    area_/=max_area;
    
    return area_;
  }

  bool ROC::target(const size_t i) const
  {
    return vec_pair_[i].first;
  }

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

}}} // of namespace statistics, yat, and theplu