source: branches/0.13-stable/test/roc.cc @ 3437

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

/*
  Copyright (C) 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 "Suite.h"

#include "yat/classifier/DataLookupWeighted1D.h"
#include "yat/random/random.h"
#include "yat/classifier/Target.h"
#include "yat/statistics/Averager.h"
#include "yat/statistics/Fisher.h"
#include "yat/statistics/ROC.h"
#include "yat/statistics/utility.h"
#include "yat/utility/Vector.h"

#include <gsl/gsl_cdf.h>

#include <cassert>
#include <cmath>
#include <deque>
#include <fstream>
#include <iostream>
#include <set>
#include <vector>

using namespace theplu::yat;

void test_empty(test::Suite&);
void test_p_approx_weighted(test::Suite& suite);
void test_p_approx_with_ties(test::Suite& suite);
void test_p_approx(test::Suite& suite);
void test_p_double_weights(test::Suite& suite);
void test_p_exact(test::Suite& suite);
void test_p_exact_weighted(test::Suite& suite);
void test_p_exact_with_ties(test::Suite& suite);
void test_p_with_weights(test::Suite& suite);
void test_p_with_weights_and_ties(test::Suite& suite);
void test_remove(test::Suite& suite);
void test_ties(test::Suite& suite);

int main(int argc, char* argv[])
{
  test::Suite suite(argc, argv);

  suite.err() << "testing ROC" << std::endl;
  utility::Vector value(31);
  std::vector<std::string> label(31,"negative");
  for (size_t i=0; i<16; i++)
    label[i] = "positive";
  classifier::Target target(label);
  for (size_t i=0; i<value.size(); i++)
    value(i)=i;
  statistics::ROC roc;
  add(roc, value.begin(), value.end(), target);
  double area = roc.area();
  if (!suite.equal(area,0.0)){
    suite.err() << "test_roc: area is " << area << " should be 0.0"
                << std::endl;
    suite.add(false);
  }
  target.set_binary(0,false);
  target.set_binary(1,true);
  roc.reset();
  add(roc, value.begin(), value.end(), target);
  area = roc.area();
  if (!suite.equal(area,1.0)){
    suite.err() << "test_roc: area is " << area << " should be 1.0"
                << std::endl;
    suite.add(false);
  }

  double p = roc.p_right();
  double p2 = roc.p_value();
  double p_matlab = 0.00000115;
  if (!(p/p_matlab < 1.01 && p/p_matlab > 0.99)){
    suite.err() << "get_p_approx: p-value not correct" << std::endl;
    suite.err() << p << " expected " << p_matlab << std::endl;
    suite.add(false);
  }
  if (!(p2==2*p)) {
    suite.add(false);
    suite.err() << "Two-sided P-value should equal 2 * p_right.\n";
  }
  roc.minimum_size() = 20;
  p = roc.p_right();
  p2 = roc.p_value();
  if (!( p < 1e-8 && p > 1e-9) ){
    suite.err() << "get_p_exact: p-value not correct" << std::endl;
    suite.add(false);
  }
  if (!( p2==2*p)) {
    suite.add(false);
    suite.err() << "Two-sided P-value should equal 2 * p_right.\n";
  }

  classifier::DataLookupWeighted1D dlw(target.size(),1.3);
  add(roc, dlw.begin(), dlw.end(), target);
  test_ties(suite);
  test_p_approx_with_ties(suite);
  test_p_exact_with_ties(suite);
  test_p_approx(suite);
  test_p_exact(suite);
  test_empty(suite);
  test_p_with_weights(suite);
  test_remove(suite);
  return suite.return_value();
}


void test_p_exact_with_ties(test::Suite& suite)
{
  suite.out() << "test p exact with ties\n";
  statistics::ROC roc;
  /*
    +++-- 6
    ++-+- 5 4.5 *** our case ***
    +-++- 4 4.5
    ++--+ 4 3.5
    +-+-+ 3 3.5
    +--++ 2 2
    -+++- 3 3
    -++-+ 2 2
    -+-++ 1 0.5 *** our second case ***
    --+++ 0 0.5
   */
  roc.add(2, true);
  roc.add(1, true);
  roc.add(1, false);
  roc.add(0, true);
  roc.add(-1, false);
  roc.area();
  if (!suite.equal(roc.p_right(), 3.0/10.0)) {
    suite.add(false);
    suite.out() << "  p_right: expected 0.3\n";
  }
  else
    suite.add(true);
  if (!suite.equal(roc.p_value(), 5.0/10.0)) {
    suite.add(false);
    suite.out() << "  (two-sided) p_value: expected 0.5\n";
  }
  else
    suite.add(true);

  suite.out() << "test p exact with ties II\n";
  roc.reset();
  roc.add(2, false);
  roc.add(1, true);
  roc.add(1, false);
  roc.add(0, true);
  roc.add(-1, true);
  suite.add(suite.equal(roc.area(), 0.5/6));
  if (!suite.add(suite.equal(roc.p_right(), 10.0/10.0)))
    suite.out() << "  p_right: expected 0.3\n";
  if (!suite.add(suite.equal(roc.p_value(), 3.0/10.0)))
    suite.out() << "  (two-sided) p_value: expected 0.5\n";
}


void test_p_approx_with_ties(test::Suite& suite)
{
  suite.out() << "test p approx with ties\n";
  statistics::ROC roc;
  for (size_t i=0; i<100; ++i) {
    roc.add(1, i<60);
    roc.add(0, i<40);
  }
  suite.add(suite.equal(roc.area(), 0.6));
  // Having only two data values, 0 and 1, data can be represented as
  // a 2x2 contigency table, and ROC test is same as Fisher's exact
  // test.
  statistics::Fisher fisher;
  fisher.oddsratio(60, 40, 40, 60);
  suite.add(suite.equal_fix(roc.p_value(), fisher.p_value(), 0.0002));
}

void test_ties(test::Suite& suite)
{
  suite.out() << "test ties\n";
  statistics::ROC roc;
  for (size_t i=0; i<20; ++i)
    roc.add(10.0, i<10);
  if (!suite.add(suite.equal(roc.area(), 0.5))) {
    suite.err() << "error: roc with ties: area: " << roc.area() << "\n";
  }
}

void test_p_exact(test::Suite& suite)
{
  suite.out() << "test_p_exact\n";
  statistics::ROC roc;
  for (size_t i=0; i<9; ++i)
    roc.add(i, i<5);
  if (roc.p_right()<0.5) {
    suite.add(false);
    suite.err() << "error: expected p-value>0.5\n  found: "
                << roc.p_right() << "\n";
  }
}


void test_p_approx(test::Suite& suite)
{
  suite.out() << "test_p_approx\n";
  statistics::ROC roc;
  for (size_t i=0; i<100; ++i)
    roc.add(i, i<50);
  if (roc.p_right()<0.5) {
    suite.add(false);
    suite.err() << "error: expected p-value>0.5\n  found: "
                << roc.p_right() << "\n";
  }
  if (roc.p_value() > 1.0) {
    suite.err() << "error: expected p-value <= 1\n    found: "
                << roc.p_value() << "\n";
    suite.add(false);
  }
}


void test_empty(test::Suite& suite)
{
  suite.err() << "test empty\n";
  // testing bug #669
  statistics::ROC roc;
  roc.p_value();
  roc.area();
  suite.err() << "test empty done\n";
}


void test_p_with_weights(test::Suite& suite)
{
  suite.out() << "test p with weights\n";
  test_p_double_weights(suite);
  test_p_exact_weighted(suite);
  test_p_approx_weighted(suite);
  test_p_with_weights_and_ties(suite);
}


void test_p_with_weights_and_ties(test::Suite& suite)
{
  suite.out() << "test p with weights and ties\n";
  statistics::ROC roc;
  roc.add(10, true, 1.0);
  roc.add(10, false, 3.0);
  roc.add(20, true, 2.0);
  roc.add(30, true, 1.0);
  if (!suite.equal(roc.area(), 0.875)) {
    suite.add(false);
    suite.out() << "roc area: " << roc.area() << "\n";
  }
  double p = roc.p_right();
  if (!suite.equal(p, 8.0/24.0)) {
    suite.add(false);
    suite.out() << "p_right() failed\n";
  }
  p = roc.p_value();
  if (!suite.equal(p, (8.0+6.0)/24.0)) {
    suite.add(false);
    suite.out() << "p_value() failed\n";
  }
}


void test_p_double_weights(test::Suite& suite)
{
  suite.out() << "test p with double weights\n";
  std::vector<double> w(5,1.0);
  w[0]=0.1;
  w[4]=10;
  std::vector<double> x(5);
  for (size_t i=0; i<x.size(); ++i)
    x[i] = i;
  statistics::ROC roc;
  statistics::ROC roc2;
  for (size_t i=0; i<x.size(); ++i) {
    roc.add(x[i], i<2, w[i]);
    roc2.add(x[i], i<2, 2*w[i]);
  }
  if (!suite.equal(roc.area(), roc2.area())) {
    suite.add(false);
    suite.err() << "area failed\n";
  }
  if (!suite.equal(roc.p_right(), roc2.p_right())) {
    suite.add(false);
    suite.err() << "exact p failed\n";
  }
  roc.minimum_size() = 0;
  roc2.minimum_size() = 0;
  if (!suite.equal(roc.p_right(), roc2.p_right())) {
    suite.add(false);
    suite.err() << "approximative p failed\n";
  }
}


void test_p_exact_weighted(test::Suite& suite)
{
  suite.out() << "test p exact weighted\n";
  std::vector<double> x(10);
  std::vector<std::pair<bool, double> > w(10);
  for (size_t i=0; i<x.size(); ++i) {
    x[i] = i;
    w[i].first = false;
    w[i].second = 1.0;
  }
  w[3].first = true;
  w[7].first = true;
  w[8].first = true;
  w[1].second = 10;
  w[7].second = 10;
  w[9].second = 0.1;

  statistics::ROC roc;
  for (size_t i=0; i<x.size(); ++i)
    roc.add(x[i], w[i].first, w[i].second);
  roc.minimum_size() = 100;

  std::sort(w.begin(), w.end());
  unsigned long perm = 0;
  unsigned long k = 0;
  unsigned long k2 = 0;
  while (true) {
    ++perm;
    statistics::ROC roc2;
    for (size_t i=0; i<x.size(); ++i)
      roc2.add(x[i], w[i].first, w[i].second);
    if (roc2.area() >= roc.area())
      ++k;
    if (roc2.area() <= 1-roc.area()+1e-10)
      ++k2;

    if (!next_permutation(w.begin(), w.end()))
      break;
  }
  double p_value = roc.p_right();
  if (!suite.add(suite.equal(p_value, static_cast<double>(k)/perm))) {
    suite.out() << "area: " << roc.area() << "\n"
                << perm << " permutations of which\n"
                << k << " with larger (or equal) area "
                << "corresponding to P=" << static_cast<double>(k)/perm << "\n"
                << "p_right() returned: " << p_value
                << "\n";
  }
  p_value = roc.p_value();
  if (!suite.add(suite.equal(p_value, static_cast<double>(k+k2)/perm))) {
    suite.out() << "area: " << roc.area() << "\n"
                << perm << " permutations of which\n"
                << k << " with larger (or equal) area and\n"
                << k2 << " with smaller (or equal) area\n"
                << "corresponding to P="
                << static_cast<double>(k+k2)/perm << "\n"
                << "p_value() returned: " << p_value
                << "\n";
  }
}


void test_p_approx_weighted(test::Suite& suite)
{
  suite.out() << "test p approx weighted\n";
  std::vector<double> x(200);
  std::vector<double> w(200, 1.0);
  std::deque<bool> label(200);

  for (size_t i=0; i<x.size(); ++i) {
    x[i] = i;
    label[i] = i>30 && i<70;
    if (i<100)
      w[i] = 100.0 / (100+i);
    else
      w[i] = 0.0001;
  }

  statistics::ROC roc;
  for (size_t i=0; i<x.size(); ++i)
    roc.add(x[i], label[i], w[i]);
  roc.minimum_size() = 0;
  double p = roc.p_right();

  std::set<size_t> checkpoints;
  size_t perm = 100000;
  checkpoints.insert(10);
  checkpoints.insert(100);
  checkpoints.insert(1000);
  checkpoints.insert(10000);
  checkpoints.insert(perm);
  statistics::Averager averager;
  for (size_t i=1; i<=perm; ++i) {
    theplu::yat::random::random_shuffle(x.begin(), x.end());
    statistics::ROC roc2;
    for (size_t j=0; j<x.size(); ++j)
      roc2.add(x[j], label[j], w[j]);
    if (roc2.area()>=roc.area())
      averager.add(1.0);
    else
      averager.add(0.0);
    if (checkpoints.find(i)!=checkpoints.end()) {
      if (gsl_cdf_binomial_P(averager.sum_x(), p, averager.n())<1e-10 ||
          gsl_cdf_binomial_Q(averager.sum_x(), p, averager.n())<1e-10) {
        suite.err() << "error: approx p value and permutation p-value "
                    << "deviate more than expected\n"
                    << "area: " << roc.area() << "\n"
                    << "approx p: " << p << "\n"
                    << "permutations: " << averager.n() << "\n"
                    << "successful: " << averager.sum_x() << "\n"
                    << "corresponds to P=" << averager.mean() << "\n";
        suite.add(false);
        return;
      }
    }
  }
}

void test_remove(test::Suite& suite)
{
  using statistics::ROC;
  ROC roc;
  roc.add(1, true);
  roc.add(2, false);
  ROC roc2(roc);
  if (!suite.add(suite.equal(roc.area(), roc2.area())))
    suite.out() << "test_remove failed: copy failed\n";
  roc.add(2.3, true, 1.2);
  try {
    roc.remove(2.3, true, 1.2);
  }
  catch (std::runtime_error& e) {
    suite.add(false);
    suite.out() << "exception what(): " << e.what() << "\n";
  }
  if (!suite.add(suite.equal(roc.area(), roc2.area())))
    suite.out() << "test remove failed\n";
  try {
    roc.remove(2, true);
    suite.out() << "no exception thrown\n";
    suite.add(false);
  }
  catch (std::runtime_error& e) {
    suite.add(true);
  }
}