source: trunk/lib/statistics/Pearson.cc @ 514

// $Id: Pearson.cc 514 2006-02-20 09:45:34Z peter $


#include <c++_tools/statistics/Pearson.h>
#include <c++_tools/statistics/AveragerPair.h>
#include <c++_tools/statistics/AveragerPairWeighted.h>
#include <c++_tools/gslapi/vector.h>
#include <c++_tools/classifier/Target.h>

#include <cmath>
#include <gsl/gsl_cdf.h>


namespace theplu {
namespace statistics {  

  Pearson::Pearson(bool b) 
    : Score(b), r_(0), nof_samples_(0)
  {
  }

  double Pearson::p_value() const
  {
    if(weighted_)
      return 1;
    if(nof_samples_<2){
      std::cerr << "Warning: Only " << nof_samples_ << "samples. " 
                << "Need at lest 3.\n";
      return 1;
    }

    double t = sqrt(nof_samples_ - 2)*fabs(r_) /sqrt(1-r_*r_);
    double p = gsl_cdf_tdist_Q(t, nof_samples_ -2 );
    if (absolute_)
      return 2*p;
    if (r_<0)
      return 1-p;
    return p;

  }

  double Pearson::score(const classifier::Target& target, 
                        const gslapi::vector& value)
  {
    weighted_=false;
    AveragerPair ap;
    for (size_t i=0; i<target.size(); i++){
      if (target.binary(i))
        ap.add(1, value(i));
      else
        ap.add(-1, value(i));
      nof_samples_ = target.size();
    }
    r_ = ap.correlation();
    if (r_<0 && absolute_)
      return -r_;
     
    return r_;
  } 
    
  double Pearson::score(const classifier::Target& target, 
                        const gslapi::vector& value,
                        const gslapi::vector& weight)
  {
    weighted_=true;
    AveragerPairWeighted ap;
    for (size_t i=0; i<target.size(); i++){
      if (target.binary(i))
        ap.add(1, value(i),1,weight(i));
      else
        ap.add(-1, value(i),1,weight(i));
      nof_samples_ = target.size();
    }
    r_ = ap.correlation();
    if (r_<0 && absolute_)
      return -r_;
     
    return r_;
  }

}} // of namespace statistics and namespace theplu