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

// $Id: ROC.cc 1797 2009-02-12 18:07:10Z peter $

/*
  Copyright (C) 2004, 2005 Peter Johansson
  Copyright (C) 2006, 2007, 2008 Jari Häkkinen, Peter Johansson

  This file is part of the yat library, http://dev.thep.lu.se/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 3 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 yat. If not, see <http://www.gnu.org/licenses/>.
*/

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

#include <gsl/gsl_cdf.h>

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

namespace theplu {
namespace yat {
namespace statistics {  

  ROC::ROC(void) 
    :minimum_size_(10)
  {
    reset();
  }

  ROC::~ROC(void)
  {
  }


  void ROC::add(double x, bool target, double w)
  {
    if (!w)
      return;
    multimap_.insert(std::make_pair(x, std::make_pair(target, w)));
    if (target)
      w_pos_+=w;
    else
      w_neg_+=w;
    area_=-1;
  }


  double ROC::area(void)
  {
    if (area_==-1){
      AUC auc(false);
      area_=auc.score(multimap_);
    }
    return area_;
  }


  double ROC::get_p_approx(double x) const
  {
    // make x standard normal
    x -= 0.5;
    // Not integrating from the middle of the bin, but from the inner edge.
    if (x>0)
      x -= 0.5/(n_pos()*n_neg());
    else if(x<0)
      x += 0.5/(n_pos()*n_neg());
    else if (x==0)
      return 0.5;
    double sigma = std::sqrt( ( n()+1.0 )/(12*n_pos()*n_neg()) );
    return gsl_cdf_gaussian_Q(x, sigma);
  }

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

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


  const unsigned int& ROC::minimum_size(void) const
  {
    return minimum_size_;
  }


  double ROC::n(void) const
  {
    return n_pos()+n_neg();
  }


  double ROC::n_neg(void) const
  {
    return w_neg_;
  }


  double ROC::n_pos(void) const
  {
    return w_pos_;
  }


  double ROC::p_value() const
  {
    double p = 2*p_value_one_sided();
    return std::min(p, 2-p); 
  }


  double ROC::p_value_one_sided() const
  {
    double area(area_);
    if (area_==-1){
      AUC auc(false);
      area = auc.score(multimap_);
    }
    if (n_pos() < minimum_size_ && n_neg() < minimum_size_) {
      // for small areas we calculate probabilitu to get larger area -
      // not larger or equal.
      if (area<0.5)
        return 1-get_p_exact((1-area)*n_pos()*n_neg(), n_pos(), n_neg());
      return get_p_exact(area*n_pos()*n_neg(), n_pos(), n_neg());
    }
    return get_p_approx(area);
  }


  void ROC::reset(void)
  {
    area_=-1;
    w_pos_=0;
    w_neg_=0;
    multimap_.clear();
  }

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