source: branches/0.13-stable/yat/statistics/ROC.cc @ 3437

// $Id: ROC.cc 3437 2015-11-20 04:44:02Z peter $

/*
  Copyright (C) 2004, 2005 Peter Johansson
  Copyright (C) 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
  Copyright (C) 2011, 2012, 2013, 2015 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 <config.h>

#include "ROC.h"
#include "AUC.h"

#include "yat/utility/Exception.h"

#include <gsl/gsl_cdf.h>

#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include <sstream>
#include <utility>

namespace theplu {
namespace yat {
namespace statistics {

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


  void ROC::add(double x, bool target, double w)
  {
    if (!w)
      return;
    ROC::Map::value_type element(x, std::make_pair(target, w));
    ROC::Map::iterator lower = multimap_.lower_bound(x);
    if (lower!=multimap_.end() && lower->first == x)
      has_ties_ = true;
    multimap_.insert(lower, element);
    if (target)
      pos_weights_.add(w);
    else
      neg_weights_.add(w);
    area_ = std::numeric_limits<double>::quiet_NaN();
  }


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


  double ROC::get_p_approx(double x) const
  {
    size_t n_pos = nof_points(pos_weights_);
    size_t n_neg = nof_points(neg_weights_);
    size_t nof_samples = n_pos + n_neg;
    // 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
      return 0.5;
    double var = 1.0+nof_samples;
    if (has_ties_) {
      double correction = 0;
      Map::const_iterator first = multimap_.begin();
      Map::const_iterator last = multimap_.begin();
      while (first!=multimap_.end()) {
        size_t n = 0;
        while (first->first == last->first) {
          ++n;
          ++last;
        }
        correction += n * (n-1) * (n+1);
        first = last;
      }
      /*
        mn(N+1)/12-[mn/(12N(N-1)) * sum(t(t-1)(t+1))] =
        mn/12 [ N+1 - 1/(N(N-1)) * sum(t(t-1)(t+1)) ]
      */
      var -= correction/(nof_samples * (nof_samples-1));
    }
    var = var / (12*n_pos*n_neg);
    return gsl_cdf_gaussian_Q(x, std::sqrt(var));
  }


  bool ROC::is_weighted(void) const
  {
    return pos_weights_.variance() || neg_weights_.variance()
      || pos_weights_.mean() != neg_weights_.mean();
  }


  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 neg_weights_.sum_x();
  }


  double ROC::n_pos(void) const
  {
    return pos_weights_.sum_x();
  }


  size_t ROC::nof_points(const Averager& a) const
  {
    return static_cast<size_t>(a.sum_x()*a.sum_x()/a.sum_xx());
  }


  double ROC::p_exact_left(double area) const
  {
    if (is_weighted())
      return p_left_weighted(area);
    return p_exact_with_ties(multimap_.rbegin(), multimap_.rend(),
                             (1-area)*pos_weights_.n()*neg_weights_.n(),
                             pos_weights_.n(), neg_weights_.n());
  }


  double ROC::p_exact_right(double area) const
  {
    if (is_weighted())
      return p_right_weighted(area);
    return p_exact_with_ties(multimap_.begin(), multimap_.end(),
                             area*pos_weights_.n()*neg_weights_.n(),
                             pos_weights_.n(), neg_weights_.n());
  }


  double ROC::p_left_weighted(double area) const
  {
    return count(utility::pair_first_iterator(multimap_.begin()),
                 utility::pair_first_iterator(multimap_.end()), 1-area);
  }


  double ROC::p_right_weighted(double area) const
  {
    return count(utility::pair_first_iterator(multimap_.rbegin()),
                 utility::pair_first_iterator(multimap_.rend()), area);
  }


  double ROC::p_left() const
  {
    if (std::isnan(area()))
      return std::numeric_limits<double>::quiet_NaN();
    if (use_exact_method())
      return p_exact_left(area());
    return get_p_approx(1-area());
  }


  double ROC::p_right() const
  {
    if (std::isnan(area()))
      return std::numeric_limits<double>::quiet_NaN();
    if (use_exact_method())
      return p_exact_right(area());
    return get_p_approx(area());
  }


  double ROC::p_value() const
  {
    if (std::isnan(area()))
      return std::numeric_limits<double>::quiet_NaN();
    if (use_exact_method()) {
      double p = 0;
      double abs_area = std::max(area(), 1-area());
      p = p_exact_right(abs_area);
      if (has_ties_) {
        p += p_exact_left(1.0 - abs_area);
      }
      else
        p *= 2.0;
      // avoid double counting when area is 0.5
      return std::min(p, 1.0);
    }
    return 2*get_p_approx(std::max(area(), 1-area()));
  }


  double ROC::p_value_one_sided() const
  {
    return p_right();
  }


  void ROC::remove(double value, bool target, double weight)
  {
    if (!weight)
      return;
    std::pair<Map::iterator, Map::iterator> iter = multimap_.equal_range(value);
    while (iter.first!=iter.second) {
      if (iter.first->second.first==target && iter.first->second.first==target){
        multimap_.erase(iter.first);
        if (target)
          pos_weights_.add(weight, -1);
        else
          neg_weights_.add(weight, -1);
        area_ = std::numeric_limits<double>::quiet_NaN();
        return;
      }
      ++iter.first;
    }
    std::stringstream ss;
    ss << "ROC::remove(" << value << "," << target << "," << weight << "): "
       << "no such element";
    throw utility::runtime_error(ss.str());
  }


  void ROC::reset(void)
  {
    area_ = std::numeric_limits<double>::quiet_NaN();
    has_ties_ = false;
    neg_weights_.reset();
    pos_weights_.reset();
    multimap_.clear();
  }


  bool ROC::use_exact_method(void) const
  {
    return (n_pos() < minimum_size_) || (n_neg() < minimum_size_);
  }


  ROC::Weights::Weights(void)
    : small_pos(0), small_neg(0), tied_pos(0), tied_neg(0)
  {}

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