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

#ifndef _theplu_yat_statistics_roc_
#define _theplu_yat_statistics_roc_

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

/*
  Copyright (C) 2004 Peter Johansson
  Copyright (C) 2005, 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 "Averager.h"
#include "yat/utility/deprecate.h"
#include "yat/utility/stl_utility.h"
#include "yat/utility/yat_assert.h"

#include <gsl/gsl_randist.h>

#include <algorithm>
#include <map>
#include <utility>
#include <vector>

namespace theplu {
namespace yat {
namespace statistics {

  ///
  /// @brief Reciever Operating Characteristic.
  ///
  /// As the area under an ROC curve is equivalent to Mann-Whitney U
  /// statistica, this class can be used to perform a Mann-Whitney
  /// U-test (aka Wilcoxon).
  ///
  /// \see AUC
  ///
  class ROC
  {

  public:
    ///
    /// @brief Default constructor
    ///
    ROC(void);

    /**
       \brief Add a data value.

       \param value data value
       \param target \c true if value belongs to class positive

       \param weight indicating how important the data point is. A
       zero weight implies the data point is ignored. A negative
       weight should be understood as removing a data point and thus
       typically only makes sense if there is a previously added data
       point with same \a value and \a target.

    */
    void add(double value, bool target, double weight=1.0);

    /**
       \brief Area Under Curve, AUC

       \see AUC for how the area is calculated

       @return Area under curve.
    */
    double area(void) const;

    /**
       \brief threshold for p_value calculation

       Function can used to change the minimum_size.

       \return reference to threshold minimum size
     */
    unsigned int& minimum_size(void);

    /**
       \brief threshold for p_value calculation

       Threshold deciding whether p-value is computed using exact
       method or a Gaussian approximation. If either number of positive
       samples, n_pos(void), or number of negative samples,
       n_neg(void), are smaller than minimum_size the exact method is
       used.

       \see p_value

       \return const reference to minimum_size
    */
    const unsigned int& minimum_size(void) const;

    ///
    /// \brief number of samples
    ///
    /// @return sum of weights
    ///
    double n(void) const;

    ///
    /// \brief number of negative samples
    ///
    /// @return sum of weights with negative target
    ///
    double n_neg(void) const;

    ///
    /// \brief number of positive samples
    ///
    /// @return sum of weights with positive target
    ///
    double n_pos(void) const;


    /**
       Calculates the probability to get this area (or less).

       \see p_right for more details
     */
    double p_left(void) const;

    /**
       \brief One-sided P-value

       Calculates the one-sided p-value, i.e., probability to get this
       area (or greater) given that there is no difference
       between the two classes.

       \b Exact \b method: In the exact method the function goes
       through all permutations and counts what fraction for which the
       area is greater (or equal) than area in original
       permutation. In case all non-zero weights are not equal,
       iterating through all permutations is not sufficient so
       algorithm goes through all combinations instead which quickly
       becomes a large number (N!).

       \b Large-sample \b Approximation: When many data points are
       available, see minimum_size(), a Gaussian approximation is used
       and the p-value is calculated as
       \f[
       P = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^z
       \exp{\left(-\frac{t^2}{2}\right)} dt
       \f]

       where

       \f[
       z = \frac{\textrm{area} - 0.5 - 0.5/(n^+ \cdot n^-)}{s}
       \f]

       and

       \f[
       s^2 = \frac{n+1+\sum \left(n_x \cdot (n_x^2-1)\right)}
       {12\cdot n^+\cdot n^-}
       \f]

       where sum runs over different data values (of ties) and \f$ n_x
       \f$ is number data points with that value. The sum is a
       correction term for ties and is zero if there are no ties.

       The number of samples in a group, \f$ n^+ \f$, is calculated as
       \f$ n = (\sum w)^2 / \sum w^2 \f$

       \return \f$ P(a \ge \textrm{area}) \f$
     */
    double p_right(void) const;

    /**
       \deprecated Provided for backward compatibility with 0.10
       API. Use p_right() instead.
     */
    double p_value_one_sided(void) const YAT_DEPRECATE;

    /**
       \brief Two-sided p-value.

       Calculates the probability to get an area, \c a, equal or more
       extreme than \c area
       \f[
       P(a \ge \textrm{max}(\textrm{area},1-\textrm{area})) +
       P(a \le \textrm{min}(\textrm{area}, 1-\textrm{area})) \f]

       If there are no ties, distribution of \a a is symmetric, so if
       area is greater than 0.5, this boils down to \f$ P = 2*P(a \ge
       \textrm{area}) = 2*P_\textrm{one-sided}\f$.

       \return two-sided p-value

       \see p_right
    */
    double p_value(void) const;

    /**
       \brief remove a data value

       A data point with identical \a value, \a target, and \a weight
       must have beed added prior calling this function; else an
       exception is thrown.

       \since New in yat 0.9
     */
    void remove(double value, bool target, double weight=1.0);

    /**
       @brief Set everything to zero
    */
    void reset(void);

  private:
    typedef std::multimap<double, std::pair<bool, double> > Map;
    typedef std::vector<std::pair<bool, Map::mapped_type> > Vector;

    // struct used in count functions
    struct Weights
    {
      Weights(void);
      double small_pos;
      double small_neg;
      double tied_pos;
      double tied_neg;
    };

    /// Implemented as in MatLab 13.1
    double get_p_approx(double) const;

    /**
       return false if all non-zero weights are equal
     */
    bool is_weighted(void) const;

    /**
       return (sum x)^2 / sum x^2
     */
    size_t nof_points(const Averager& a) const;

    /*
      Calculate probability to get an area equal (smaller) than \a
      area given the distribution of weights and ties in multimap_
     */
    double p_left_weighted(double area) const;

    /*
      Calculate probability to get an area equal (greater) than \a
      area given the distribution of weights and ties in multimap_
     */
    double p_right_weighted(double area) const;

    /*
      Count number of combinations (of N!) that gives weight sum equal
      or larger than \a threshold.

      Range [first, last) is used to check for ties. If, e.g., *first
      and *(first+1) are equal implies that the two largest values are
      equal.
     */
    template <typename Iterator>
    double count(Iterator first, Iterator last, double threshold) const;

    /*
      Loop over all elements in \a weights and call count(7)
     */
    template <typename Iterator>
    double count(Vector& weights, Iterator iter, Iterator last,
                 double threshold, double sum, const Weights& weight) const;

    /*
      Count number of combinations in which sum>=threshold given
      classes and weights in \a weight. Range [iter, last) is used to
      handle ties.
     */
    template <typename Iterator>
    double count(Vector& weights, Iterator iter, Iterator last,
                 double threshold, double sum, Weights weight,
                 const std::pair<bool, double>& entry) const;

    /*
      Calculates probability to get \a block number of pairs correctly
      sorted when having \a pos positive samples and \a neg negative
      samples given the distribution of ties as in [first, last).
     */
    template<typename ForwardIterator>
    double p_exact_with_ties(ForwardIterator first, ForwardIterator last,
                             double block, unsigned int pos,
                             unsigned int neg) const;

    /**
       \return P(auc >= area)
     */
    double p_exact_right(double area) const;

    /**
       \return P(auc <= area)
     */
    double p_exact_left(double area) const;

    bool use_exact_method(void) const;

    mutable double area_;
    bool has_ties_;
    unsigned int minimum_size_;
    Averager neg_weights_;
    Averager pos_weights_;
    Map multimap_;
  };

  template<typename ForwardIterator>
  double
  ROC::p_exact_with_ties(ForwardIterator begin, ForwardIterator end,
                         double block, unsigned int pos,unsigned int neg) const
  {
    if (block <= 0)
      return 1.0;
    if (block > pos*neg)
      return 0.0;

    ForwardIterator iter(begin);
    unsigned int n=0;
    while (iter!=end && iter->first == begin->first) {
      ++iter;
      ++n;
    }
    double result = 0;
    /*
      pos1  neg1  |  n
      pos2  neg2  |
      ----  ----   ----
      pos   neg
     */

    // ensure pos1 and neg2 are non-negative
    unsigned int pos1 = n - std::min(n, neg);
    // ensure pos2 and neg1 are non-negative
    unsigned int max = std::min(n, pos);
    YAT_ASSERT(pos1<=max);
    for ( ; pos1<=max; ++pos1) {
      unsigned int neg1 = n-pos1;
      YAT_ASSERT(neg1<=n);
      unsigned int pos2 = pos-pos1;
      YAT_ASSERT(pos2<=pos);
      unsigned int neg2 = neg-neg1;
      YAT_ASSERT(neg2<=neg);
      result += gsl_ran_hypergeometric_pdf(pos1, static_cast<unsigned int>(pos),
                                           static_cast<unsigned int>(neg), n)
        * p_exact_with_ties(iter, end,
                            block - pos2*neg1 - 0.5*pos1*neg1,
                            pos2, neg2);
    }
    return result;
  }


  template <typename Iterator>
  double ROC::count(Iterator first, Iterator last, double threshold) const
  {
    Vector vec;
    vec.reserve(multimap_.size());
    // copy values from multimap_ to v
    for (Map::const_iterator i = multimap_.begin(); i!=multimap_.end(); ++i)
      vec.push_back(std::make_pair(false, i->second));

    ROC::Weights w;
    w.small_pos = pos_weights_.sum_x();
    w.small_neg = neg_weights_.sum_x();
    return count(vec, first, last, threshold*w.small_pos*w.small_neg, 0, w);
  }



  template <typename Iterator>
  double ROC::count(ROC::Vector& v, Iterator iter, Iterator last,
                    double threshold, double sum, const Weights& w) const
  {
    double result = 0.0;
    // loop over all elements
    int nof_elements = 0;
    for (ROC::Vector::iterator i=v.begin(); i!=v.end(); ++i) {
      if (i->first)
        continue;
      i->first = true;
      result += count(v, iter, last, threshold, sum, w, i->second);
      i->first = false;
      ++nof_elements;
    }
    YAT_ASSERT(nof_elements);
    return result/nof_elements;
  }


  template <typename Iterator>
  double ROC::count(Vector& weights, Iterator iter, Iterator last,
                    double threshold, double sum, Weights w,
                    const std::pair<bool, double>& entry) const
  {
    double tiny = 10e-10;

    Iterator next(iter);
    YAT_ASSERT(next!=last);
    ++next;

    // update weights
    if (entry.first) {
      w.tied_pos += entry.second;
      w.small_pos -= entry.second;
    }
    else {
      w.tied_neg += entry.second;
      w.small_neg -= entry.second;
    }

    // last entry in equal range
    if (next==last || *next!=*iter) {
      sum += 0.5*w.tied_pos*w.tied_neg + w.tied_pos * w.small_neg;
      w.tied_pos=0;
      w.tied_neg=0;
    }

    // max sum happens if all pos values belong to current equal range
    // and none of the remaining neg values
    double max_sum = sum + 0.5*(w.tied_pos+w.small_pos)*w.tied_neg +
      (w.tied_pos+w.small_pos)*w.small_neg;

    if (max_sum<threshold-tiny)
      return 0.0;
    if (sum + 0.5*w.tied_pos*(w.tied_neg+w.small_neg) >= threshold-tiny)
      return 1.0;

    if (next!=last)
      return count(weights, next, last, threshold, sum, w);
    return 0.0;
  }

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