yat/statistics/ROC.h

Code
Comments
Other
Rev Date Author Line
2696 28 Feb 12 peter 1 #ifndef _theplu_yat_statistics_roc_
2696 28 Feb 12 peter 2 #define _theplu_yat_statistics_roc_
69 29 Apr 04 peter 3
616 31 Aug 06 jari 4 // $Id$
616 31 Aug 06 jari 5
675 10 Oct 06 jari 6 /*
831 27 Mar 07 peter 7   Copyright (C) 2004 Peter Johansson
2119 12 Dec 09 peter 8   Copyright (C) 2005, 2006, 2007, 2008 Jari Häkkinen, Peter Johansson
3439 20 Nov 15 peter 9   Copyright (C) 2011, 2012, 2013, 2015 Peter Johansson
69 29 Apr 04 peter 10
1437 25 Aug 08 peter 11   This file is part of the yat library, http://dev.thep.lu.se/yat
675 10 Oct 06 jari 12
675 10 Oct 06 jari 13   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 14   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 15   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 16   License, or (at your option) any later version.
675 10 Oct 06 jari 17
675 10 Oct 06 jari 18   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 19   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 20   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 21   General Public License for more details.
675 10 Oct 06 jari 22
675 10 Oct 06 jari 23   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 24   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 25 */
675 10 Oct 06 jari 26
2709 15 Mar 12 peter 27 #include "Averager.h"
3007 01 Apr 13 peter 28 #include "yat/utility/deprecate.h"
2718 12 Apr 12 peter 29 #include "yat/utility/stl_utility.h"
2710 16 Mar 12 peter 30 #include "yat/utility/yat_assert.h"
2709 15 Mar 12 peter 31
2592 27 Oct 11 peter 32 #include <gsl/gsl_randist.h>
2592 27 Oct 11 peter 33
2736 29 May 12 peter 34 #include <algorithm>
821 18 Mar 07 peter 35 #include <map>
69 29 Apr 04 peter 36 #include <utility>
2736 29 May 12 peter 37 #include <vector>
69 29 Apr 04 peter 38
69 29 Apr 04 peter 39 namespace theplu {
680 11 Oct 06 jari 40 namespace yat {
2696 28 Feb 12 peter 41 namespace statistics {
414 01 Dec 05 peter 42
69 29 Apr 04 peter 43   ///
2696 28 Feb 12 peter 44   /// @brief Reciever Operating Characteristic.
2594 30 Oct 11 peter 45   ///
589 24 Aug 06 peter 46   /// As the area under an ROC curve is equivalent to Mann-Whitney U
589 24 Aug 06 peter 47   /// statistica, this class can be used to perform a Mann-Whitney
589 24 Aug 06 peter 48   /// U-test (aka Wilcoxon).
589 24 Aug 06 peter 49   ///
2594 30 Oct 11 peter 50   /// \see AUC
2594 30 Oct 11 peter 51   ///
779 05 Mar 07 peter 52   class ROC
69 29 Apr 04 peter 53   {
2696 28 Feb 12 peter 54
69 29 Apr 04 peter 55   public:
69 29 Apr 04 peter 56     ///
703 18 Dec 06 jari 57     /// @brief Default constructor
102 15 Jun 04 peter 58     ///
779 05 Mar 07 peter 59     ROC(void);
2696 28 Feb 12 peter 60
821 18 Mar 07 peter 61     /**
2594 30 Oct 11 peter 62        \brief Add a data value.
2594 30 Oct 11 peter 63
2594 30 Oct 11 peter 64        \param value data value
2594 30 Oct 11 peter 65        \param target \c true if value belongs to class positive
2594 30 Oct 11 peter 66
2594 30 Oct 11 peter 67        \param weight indicating how important the data point is. A
2594 30 Oct 11 peter 68        zero weight implies the data point is ignored. A negative
2594 30 Oct 11 peter 69        weight should be understood as removing a data point and thus
2594 30 Oct 11 peter 70        typically only makes sense if there is a previously added data
2594 30 Oct 11 peter 71        point with same \a value and \a target.
2594 30 Oct 11 peter 72
821 18 Mar 07 peter 73     */
821 18 Mar 07 peter 74     void add(double value, bool target, double weight=1.0);
821 18 Mar 07 peter 75
821 18 Mar 07 peter 76     /**
2594 30 Oct 11 peter 77        \brief Area Under Curve, AUC
821 18 Mar 07 peter 78
2594 30 Oct 11 peter 79        \see AUC for how the area is calculated
2594 30 Oct 11 peter 80
821 18 Mar 07 peter 81        @return Area under curve.
821 18 Mar 07 peter 82     */
3023 06 Apr 13 peter 83     double area(void) const;
821 18 Mar 07 peter 84
2594 30 Oct 11 peter 85     /**
2594 30 Oct 11 peter 86        \brief threshold for p_value calculation
2594 30 Oct 11 peter 87
2594 30 Oct 11 peter 88        Function can used to change the minimum_size.
2594 30 Oct 11 peter 89
2594 30 Oct 11 peter 90        \return reference to threshold minimum size
2594 30 Oct 11 peter 91      */
1271 09 Apr 08 peter 92     unsigned int& minimum_size(void);
718 26 Dec 06 jari 93
821 18 Mar 07 peter 94     /**
2594 30 Oct 11 peter 95        \brief threshold for p_value calculation
2594 30 Oct 11 peter 96
2594 30 Oct 11 peter 97        Threshold deciding whether p-value is computed using exact
3439 20 Nov 15 peter 98        method or a Gaussian approximation. If either number of positive
3439 20 Nov 15 peter 99        samples, n_pos(void), or number of negative samples,
2594 30 Oct 11 peter 100        n_neg(void), are smaller than minimum_size the exact method is
2594 30 Oct 11 peter 101        used.
2696 28 Feb 12 peter 102
2594 30 Oct 11 peter 103        \see p_value
2594 30 Oct 11 peter 104
2594 30 Oct 11 peter 105        \return const reference to minimum_size
821 18 Mar 07 peter 106     */
1271 09 Apr 08 peter 107     const unsigned int& minimum_size(void) const;
821 18 Mar 07 peter 108
718 26 Dec 06 jari 109     ///
2594 30 Oct 11 peter 110     /// \brief number of samples
2594 30 Oct 11 peter 111     ///
821 18 Mar 07 peter 112     /// @return sum of weights
718 26 Dec 06 jari 113     ///
848 28 Apr 07 jari 114     double n(void) const;
718 26 Dec 06 jari 115
718 26 Dec 06 jari 116     ///
2594 30 Oct 11 peter 117     /// \brief number of negative samples
2594 30 Oct 11 peter 118     ///
821 18 Mar 07 peter 119     /// @return sum of weights with negative target
718 26 Dec 06 jari 120     ///
848 28 Apr 07 jari 121     double n_neg(void) const;
821 18 Mar 07 peter 122
821 18 Mar 07 peter 123     ///
2594 30 Oct 11 peter 124     /// \brief number of positive samples
2594 30 Oct 11 peter 125     ///
821 18 Mar 07 peter 126     /// @return sum of weights with positive target
821 18 Mar 07 peter 127     ///
848 28 Apr 07 jari 128     double n_pos(void) const;
718 26 Dec 06 jari 129
3006 01 Apr 13 peter 130
2594 30 Oct 11 peter 131     /**
3006 01 Apr 13 peter 132        Calculates the probability to get this area (or less).
3006 01 Apr 13 peter 133
3006 01 Apr 13 peter 134        \see p_right for more details
3006 01 Apr 13 peter 135      */
3006 01 Apr 13 peter 136     double p_left(void) const;
3006 01 Apr 13 peter 137
3006 01 Apr 13 peter 138     /**
2594 30 Oct 11 peter 139        \brief One-sided P-value
2594 30 Oct 11 peter 140
2594 30 Oct 11 peter 141        Calculates the one-sided p-value, i.e., probability to get this
2594 30 Oct 11 peter 142        area (or greater) given that there is no difference
2594 30 Oct 11 peter 143        between the two classes.
2594 30 Oct 11 peter 144
2594 30 Oct 11 peter 145        \b Exact \b method: In the exact method the function goes
2594 30 Oct 11 peter 146        through all permutations and counts what fraction for which the
2718 12 Apr 12 peter 147        area is greater (or equal) than area in original
2718 12 Apr 12 peter 148        permutation. In case all non-zero weights are not equal,
2718 12 Apr 12 peter 149        iterating through all permutations is not sufficient so
2718 12 Apr 12 peter 150        algorithm goes through all combinations instead which quickly
2718 12 Apr 12 peter 151        becomes a large number (N!).
2594 30 Oct 11 peter 152
2594 30 Oct 11 peter 153        \b Large-sample \b Approximation: When many data points are
2594 30 Oct 11 peter 154        available, see minimum_size(), a Gaussian approximation is used
2594 30 Oct 11 peter 155        and the p-value is calculated as
2594 30 Oct 11 peter 156        \f[
2696 28 Feb 12 peter 157        P = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^z
2594 30 Oct 11 peter 158        \exp{\left(-\frac{t^2}{2}\right)} dt
2594 30 Oct 11 peter 159        \f]
2594 30 Oct 11 peter 160
2696 28 Feb 12 peter 161        where
2594 30 Oct 11 peter 162
2594 30 Oct 11 peter 163        \f[
2696 28 Feb 12 peter 164        z = \frac{\textrm{area} - 0.5 - 0.5/(n^+ \cdot n^-)}{s}
2594 30 Oct 11 peter 165        \f]
2594 30 Oct 11 peter 166
2696 28 Feb 12 peter 167        and
2594 30 Oct 11 peter 168
2594 30 Oct 11 peter 169        \f[
2594 30 Oct 11 peter 170        s^2 = \frac{n+1+\sum \left(n_x \cdot (n_x^2-1)\right)}
2594 30 Oct 11 peter 171        {12\cdot n^+\cdot n^-}
2594 30 Oct 11 peter 172        \f]
2594 30 Oct 11 peter 173
2594 30 Oct 11 peter 174        where sum runs over different data values (of ties) and \f$ n_x
2709 15 Mar 12 peter 175        \f$ is number data points with that value. The sum is a
2594 30 Oct 11 peter 176        correction term for ties and is zero if there are no ties.
2594 30 Oct 11 peter 177
2709 15 Mar 12 peter 178        The number of samples in a group, \f$ n^+ \f$, is calculated as
2709 15 Mar 12 peter 179        \f$ n = (\sum w)^2 / \sum w^2 \f$
2709 15 Mar 12 peter 180
2594 30 Oct 11 peter 181        \return \f$ P(a \ge \textrm{area}) \f$
2594 30 Oct 11 peter 182      */
3006 01 Apr 13 peter 183     double p_right(void) const;
3006 01 Apr 13 peter 184
3006 01 Apr 13 peter 185     /**
3006 01 Apr 13 peter 186        \deprecated Provided for backward compatibility with 0.10
3007 01 Apr 13 peter 187        API. Use p_right() instead.
3006 01 Apr 13 peter 188      */
3007 01 Apr 13 peter 189     double p_value_one_sided(void) const YAT_DEPRECATE;
2696 28 Feb 12 peter 190
821 18 Mar 07 peter 191     /**
2594 30 Oct 11 peter 192        \brief Two-sided p-value.
821 18 Mar 07 peter 193
2594 30 Oct 11 peter 194        Calculates the probability to get an area, \c a, equal or more
2594 30 Oct 11 peter 195        extreme than \c area
2696 28 Feb 12 peter 196        \f[
2696 28 Feb 12 peter 197        P(a \ge \textrm{max}(\textrm{area},1-\textrm{area})) +
2594 30 Oct 11 peter 198        P(a \le \textrm{min}(\textrm{area}, 1-\textrm{area})) \f]
2594 30 Oct 11 peter 199
2594 30 Oct 11 peter 200        If there are no ties, distribution of \a a is symmetric, so if
2594 30 Oct 11 peter 201        area is greater than 0.5, this boils down to \f$ P = 2*P(a \ge
2594 30 Oct 11 peter 202        \textrm{area}) = 2*P_\textrm{one-sided}\f$.
2594 30 Oct 11 peter 203
2594 30 Oct 11 peter 204        \return two-sided p-value
2594 30 Oct 11 peter 205
3006 01 Apr 13 peter 206        \see p_right
650 15 Sep 06 peter 207     */
821 18 Mar 07 peter 208     double p_value(void) const;
623 05 Sep 06 peter 209
821 18 Mar 07 peter 210     /**
2720 12 Apr 12 peter 211        \brief remove a data value
2720 12 Apr 12 peter 212
2720 12 Apr 12 peter 213        A data point with identical \a value, \a target, and \a weight
2720 12 Apr 12 peter 214        must have beed added prior calling this function; else an
2720 12 Apr 12 peter 215        exception is thrown.
2720 12 Apr 12 peter 216
2720 12 Apr 12 peter 217        \since New in yat 0.9
2720 12 Apr 12 peter 218      */
2720 12 Apr 12 peter 219     void remove(double value, bool target, double weight=1.0);
2720 12 Apr 12 peter 220
2720 12 Apr 12 peter 221     /**
821 18 Mar 07 peter 222        @brief Set everything to zero
650 15 Sep 06 peter 223     */
821 18 Mar 07 peter 224     void reset(void);
179 04 Oct 04 peter 225
475 22 Dec 05 peter 226   private:
2585 25 Oct 11 peter 227     typedef std::multimap<double, std::pair<bool, double> > Map;
2736 29 May 12 peter 228     typedef std::vector<std::pair<bool, Map::mapped_type> > Vector;
2696 28 Feb 12 peter 229
2732 08 May 12 peter 230     // struct used in count functions
2718 12 Apr 12 peter 231     struct Weights
2718 12 Apr 12 peter 232     {
2718 12 Apr 12 peter 233       Weights(void);
2718 12 Apr 12 peter 234       double small_pos;
2718 12 Apr 12 peter 235       double small_neg;
2718 12 Apr 12 peter 236       double tied_pos;
2718 12 Apr 12 peter 237       double tied_neg;
2718 12 Apr 12 peter 238     };
2718 12 Apr 12 peter 239
78 04 May 04 peter 240     /// Implemented as in MatLab 13.1
821 18 Mar 07 peter 241     double get_p_approx(double) const;
295 29 Apr 05 peter 242
2709 15 Mar 12 peter 243     /**
2718 12 Apr 12 peter 244        return false if all non-zero weights are equal
2718 12 Apr 12 peter 245      */
2718 12 Apr 12 peter 246     bool is_weighted(void) const;
2718 12 Apr 12 peter 247
2718 12 Apr 12 peter 248     /**
2709 15 Mar 12 peter 249        return (sum x)^2 / sum x^2
2709 15 Mar 12 peter 250      */
2709 15 Mar 12 peter 251     size_t nof_points(const Averager& a) const;
2709 15 Mar 12 peter 252
2585 25 Oct 11 peter 253     /*
2718 12 Apr 12 peter 254       Calculate probability to get an area equal (smaller) than \a
2718 12 Apr 12 peter 255       area given the distribution of weights and ties in multimap_
2585 25 Oct 11 peter 256      */
2718 12 Apr 12 peter 257     double p_left_weighted(double area) const;
2718 12 Apr 12 peter 258
2718 12 Apr 12 peter 259     /*
2718 12 Apr 12 peter 260       Calculate probability to get an area equal (greater) than \a
2718 12 Apr 12 peter 261       area given the distribution of weights and ties in multimap_
2718 12 Apr 12 peter 262      */
2718 12 Apr 12 peter 263     double p_right_weighted(double area) const;
2718 12 Apr 12 peter 264
2718 12 Apr 12 peter 265     /*
2718 12 Apr 12 peter 266       Count number of combinations (of N!) that gives weight sum equal
2718 12 Apr 12 peter 267       or larger than \a threshold.
2718 12 Apr 12 peter 268
2718 12 Apr 12 peter 269       Range [first, last) is used to check for ties. If, e.g., *first
2718 12 Apr 12 peter 270       and *(first+1) are equal implies that the two largest values are
2718 12 Apr 12 peter 271       equal.
2718 12 Apr 12 peter 272      */
2718 12 Apr 12 peter 273     template <typename Iterator>
2718 12 Apr 12 peter 274     double count(Iterator first, Iterator last, double threshold) const;
2718 12 Apr 12 peter 275
2718 12 Apr 12 peter 276     /*
2718 12 Apr 12 peter 277       Loop over all elements in \a weights and call count(7)
2718 12 Apr 12 peter 278      */
2718 12 Apr 12 peter 279     template <typename Iterator>
2736 29 May 12 peter 280     double count(Vector& weights, Iterator iter, Iterator last,
2718 12 Apr 12 peter 281                  double threshold, double sum, const Weights& weight) const;
2718 12 Apr 12 peter 282
2718 12 Apr 12 peter 283     /*
2718 12 Apr 12 peter 284       Count number of combinations in which sum>=threshold given
2718 12 Apr 12 peter 285       classes and weights in \a weight. Range [iter, last) is used to
2718 12 Apr 12 peter 286       handle ties.
2718 12 Apr 12 peter 287      */
2718 12 Apr 12 peter 288     template <typename Iterator>
2736 29 May 12 peter 289     double count(Vector& weights, Iterator iter, Iterator last,
2718 12 Apr 12 peter 290                  double threshold, double sum, Weights weight,
2718 12 Apr 12 peter 291                  const std::pair<bool, double>& entry) const;
2718 12 Apr 12 peter 292
2718 12 Apr 12 peter 293     /*
2718 12 Apr 12 peter 294       Calculates probability to get \a block number of pairs correctly
2718 12 Apr 12 peter 295       sorted when having \a pos positive samples and \a neg negative
2718 12 Apr 12 peter 296       samples given the distribution of ties as in [first, last).
2718 12 Apr 12 peter 297      */
2585 25 Oct 11 peter 298     template<typename ForwardIterator>
2585 25 Oct 11 peter 299     double p_exact_with_ties(ForwardIterator first, ForwardIterator last,
2710 16 Mar 12 peter 300                              double block, unsigned int pos,
2710 16 Mar 12 peter 301                              unsigned int neg) const;
295 29 Apr 05 peter 302
2585 25 Oct 11 peter 303     /**
2718 12 Apr 12 peter 304        \return P(auc >= area)
2718 12 Apr 12 peter 305      */
2718 12 Apr 12 peter 306     double p_exact_right(double area) const;
2585 25 Oct 11 peter 307
2718 12 Apr 12 peter 308     /**
2718 12 Apr 12 peter 309        \return P(auc <= area)
2718 12 Apr 12 peter 310      */
2718 12 Apr 12 peter 311     double p_exact_left(double area) const;
2585 25 Oct 11 peter 312
2585 25 Oct 11 peter 313     bool use_exact_method(void) const;
2585 25 Oct 11 peter 314
3023 06 Apr 13 peter 315     mutable double area_;
2556 20 Aug 11 peter 316     bool has_ties_;
1271 09 Apr 08 peter 317     unsigned int minimum_size_;
2709 15 Mar 12 peter 318     Averager neg_weights_;
2709 15 Mar 12 peter 319     Averager pos_weights_;
2556 20 Aug 11 peter 320     Map multimap_;
69 29 Apr 04 peter 321    };
69 29 Apr 04 peter 322
2585 25 Oct 11 peter 323   template<typename ForwardIterator>
2696 28 Feb 12 peter 324   double
2585 25 Oct 11 peter 325   ROC::p_exact_with_ties(ForwardIterator begin, ForwardIterator end,
2710 16 Mar 12 peter 326                          double block, unsigned int pos,unsigned int neg) const
2585 25 Oct 11 peter 327   {
2585 25 Oct 11 peter 328     if (block <= 0)
2585 25 Oct 11 peter 329       return 1.0;
2585 25 Oct 11 peter 330     if (block > pos*neg)
2585 25 Oct 11 peter 331       return 0.0;
2585 25 Oct 11 peter 332
2585 25 Oct 11 peter 333     ForwardIterator iter(begin);
2710 16 Mar 12 peter 334     unsigned int n=0;
2585 25 Oct 11 peter 335     while (iter!=end && iter->first == begin->first) {
2585 25 Oct 11 peter 336       ++iter;
2585 25 Oct 11 peter 337       ++n;
2585 25 Oct 11 peter 338     }
2585 25 Oct 11 peter 339     double result = 0;
2710 16 Mar 12 peter 340     /*
2710 16 Mar 12 peter 341       pos1  neg1  |  n
2710 16 Mar 12 peter 342       pos2  neg2  |
2710 16 Mar 12 peter 343       ----  ----   ----
2710 16 Mar 12 peter 344       pos   neg
2710 16 Mar 12 peter 345      */
2710 16 Mar 12 peter 346
2710 16 Mar 12 peter 347     // ensure pos1 and neg2 are non-negative
2710 16 Mar 12 peter 348     unsigned int pos1 = n - std::min(n, neg);
2710 16 Mar 12 peter 349     // ensure pos2 and neg1 are non-negative
2710 16 Mar 12 peter 350     unsigned int max = std::min(n, pos);
2710 16 Mar 12 peter 351     YAT_ASSERT(pos1<=max);
2710 16 Mar 12 peter 352     for ( ; pos1<=max; ++pos1) {
2710 16 Mar 12 peter 353       unsigned int neg1 = n-pos1;
2710 16 Mar 12 peter 354       YAT_ASSERT(neg1<=n);
2710 16 Mar 12 peter 355       unsigned int pos2 = pos-pos1;
2710 16 Mar 12 peter 356       YAT_ASSERT(pos2<=pos);
2710 16 Mar 12 peter 357       unsigned int neg2 = neg-neg1;
2710 16 Mar 12 peter 358       YAT_ASSERT(neg2<=neg);
2710 16 Mar 12 peter 359       result += gsl_ran_hypergeometric_pdf(pos1, static_cast<unsigned int>(pos),
2673 03 Dec 11 peter 360                                            static_cast<unsigned int>(neg), n)
2696 28 Feb 12 peter 361         * p_exact_with_ties(iter, end,
2696 28 Feb 12 peter 362                             block - pos2*neg1 - 0.5*pos1*neg1,
2585 25 Oct 11 peter 363                             pos2, neg2);
2585 25 Oct 11 peter 364     }
2585 25 Oct 11 peter 365     return result;
2585 25 Oct 11 peter 366   }
2585 25 Oct 11 peter 367
2718 12 Apr 12 peter 368
2718 12 Apr 12 peter 369   template <typename Iterator>
2718 12 Apr 12 peter 370   double ROC::count(Iterator first, Iterator last, double threshold) const
2718 12 Apr 12 peter 371   {
2736 29 May 12 peter 372     Vector vec;
2736 29 May 12 peter 373     vec.reserve(multimap_.size());
2736 29 May 12 peter 374     // copy values from multimap_ to v
2736 29 May 12 peter 375     for (Map::const_iterator i = multimap_.begin(); i!=multimap_.end(); ++i)
2736 29 May 12 peter 376       vec.push_back(std::make_pair(false, i->second));
2736 29 May 12 peter 377
2718 12 Apr 12 peter 378     ROC::Weights w;
2718 12 Apr 12 peter 379     w.small_pos = pos_weights_.sum_x();
2718 12 Apr 12 peter 380     w.small_neg = neg_weights_.sum_x();
2736 29 May 12 peter 381     return count(vec, first, last, threshold*w.small_pos*w.small_neg, 0, w);
2718 12 Apr 12 peter 382   }
2718 12 Apr 12 peter 383
2718 12 Apr 12 peter 384
2718 12 Apr 12 peter 385
2718 12 Apr 12 peter 386   template <typename Iterator>
2736 29 May 12 peter 387   double ROC::count(ROC::Vector& v, Iterator iter, Iterator last,
2718 12 Apr 12 peter 388                     double threshold, double sum, const Weights& w) const
2718 12 Apr 12 peter 389   {
2718 12 Apr 12 peter 390     double result = 0.0;
2718 12 Apr 12 peter 391     // loop over all elements
2736 29 May 12 peter 392     int nof_elements = 0;
2736 29 May 12 peter 393     for (ROC::Vector::iterator i=v.begin(); i!=v.end(); ++i) {
2736 29 May 12 peter 394       if (i->first)
2736 29 May 12 peter 395         continue;
2736 29 May 12 peter 396       i->first = true;
2736 29 May 12 peter 397       result += count(v, iter, last, threshold, sum, w, i->second);
2736 29 May 12 peter 398       i->first = false;
2736 29 May 12 peter 399       ++nof_elements;
2718 12 Apr 12 peter 400     }
2736 29 May 12 peter 401     YAT_ASSERT(nof_elements);
2736 29 May 12 peter 402     return result/nof_elements;
2718 12 Apr 12 peter 403   }
2718 12 Apr 12 peter 404
2736 29 May 12 peter 405
2718 12 Apr 12 peter 406   template <typename Iterator>
2736 29 May 12 peter 407   double ROC::count(Vector& weights, Iterator iter, Iterator last,
2718 12 Apr 12 peter 408                     double threshold, double sum, Weights w,
2718 12 Apr 12 peter 409                     const std::pair<bool, double>& entry) const
2718 12 Apr 12 peter 410   {
2718 12 Apr 12 peter 411     double tiny = 10e-10;
2718 12 Apr 12 peter 412
2718 12 Apr 12 peter 413     Iterator next(iter);
2736 29 May 12 peter 414     YAT_ASSERT(next!=last);
2718 12 Apr 12 peter 415     ++next;
2718 12 Apr 12 peter 416
2718 12 Apr 12 peter 417     // update weights
2718 12 Apr 12 peter 418     if (entry.first) {
2718 12 Apr 12 peter 419       w.tied_pos += entry.second;
2718 12 Apr 12 peter 420       w.small_pos -= entry.second;
2718 12 Apr 12 peter 421     }
2718 12 Apr 12 peter 422     else {
2718 12 Apr 12 peter 423       w.tied_neg += entry.second;
2718 12 Apr 12 peter 424       w.small_neg -= entry.second;
2718 12 Apr 12 peter 425     }
2718 12 Apr 12 peter 426
2718 12 Apr 12 peter 427     // last entry in equal range
2718 12 Apr 12 peter 428     if (next==last || *next!=*iter) {
2718 12 Apr 12 peter 429       sum += 0.5*w.tied_pos*w.tied_neg + w.tied_pos * w.small_neg;
2718 12 Apr 12 peter 430       w.tied_pos=0;
2718 12 Apr 12 peter 431       w.tied_neg=0;
2718 12 Apr 12 peter 432     }
2718 12 Apr 12 peter 433
2718 12 Apr 12 peter 434     // max sum happens if all pos values belong to current equal range
2718 12 Apr 12 peter 435     // and none of the remaining neg values
2718 12 Apr 12 peter 436     double max_sum = sum + 0.5*(w.tied_pos+w.small_pos)*w.tied_neg +
2718 12 Apr 12 peter 437       (w.tied_pos+w.small_pos)*w.small_neg;
2718 12 Apr 12 peter 438
2718 12 Apr 12 peter 439     if (max_sum<threshold-tiny)
2718 12 Apr 12 peter 440       return 0.0;
2719 12 Apr 12 peter 441     if (sum + 0.5*w.tied_pos*(w.tied_neg+w.small_neg) >= threshold-tiny)
2718 12 Apr 12 peter 442       return 1.0;
2718 12 Apr 12 peter 443
2718 12 Apr 12 peter 444     if (next!=last)
2718 12 Apr 12 peter 445       return count(weights, next, last, threshold, sum, w);
2718 12 Apr 12 peter 446     return 0.0;
2718 12 Apr 12 peter 447   }
2718 12 Apr 12 peter 448
683 11 Oct 06 jari 449 }}} // of namespace statistics, yat, and theplu
69 29 Apr 04 peter 450 #endif

Generated on Sat Oct 14 02:23:57 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)