Changeset 623 for trunk

Timestamp:

Sep 5, 2006, 4:13:12 AM (17 years ago)

Author:

Peter

Message:

fixes #112 and refs #123 added overloaded function score taking Target and DataLookupWeighted1D, which is needed for InputRanker?.

Location:

trunk/c++_tools/statistics

Files:

• Fisher.cc (modified) (4 diffs)
• Fisher.h (modified) (1 diff)
• FoldChange.cc (modified) (2 diffs)
• FoldChange.h (modified) (2 diffs)
• Pearson.cc (modified) (2 diffs)
• Pearson.h (modified) (1 diff)
• ROC.cc (modified) (2 diffs)
• ROC.h (modified) (1 diff)
• SNR.cc (modified) (2 diffs)
• SNR.h (modified) (2 diffs)
• Score.h (modified) (2 diffs)
• WilcoxonFoldChange.cc (modified) (1 diff)
• WilcoxonFoldChange.h (modified) (1 diff)
• tScore.cc (modified) (2 diffs)
• tScore.h (modified) (1 diff)

trunk/c++_tools/statistics/Fisher.cc

@@ -4 +4 @@
 #include <c++_tools/statistics/Score.h>
 #include <c++_tools/statistics/utility.h>
+#include <c++_tools/classifier/DataLookupWeighted1D.h>
 #include <c++_tools/classifier/Target.h>
 
@@ -136 +137 @@
   }
   
-    double Fisher::score(const classifier::Target& target, 
-                         const utility::vector& value, 
-                         const utility::vector& weight) 
+  double Fisher::score(const classifier::Target& target, 
+                       const classifier::DataLookupWeighted1D& value) 
   {
     weighted_=true;
@@ -144 +144 @@
     for (size_t i=0; i<target.size(); i++)
       if (target.binary(i))
+        if (value.data(i)>value_cutoff_)
+          a_+=value.weight(i);
+        else
+          c_+=value.weight(i);
+      else
+        if (value.data(i)>value_cutoff_)
+          b_+=value.weight(i);
+        else
+          d_+=value.weight(i);
+        
+    // If a column sum or a row sum is zero, the table is non-sense
+    if ((a_==0 || d_==0) && (c_==0 || b_==0)){
+      // Peter should throw an exception here
+      std::cerr << "Warning: Fisher: Table is not valid\n";
+      return 1;
+    }
+
+    return oddsratio(a_,b_,c_,d_); 
+  }
+  
+    double Fisher::score(const classifier::Target& target, 
+                         const utility::vector& value, 
+                         const utility::vector& weight) 
+  {
+    weighted_=true;
+    a_=b_=c_=d_=0;
+    for (size_t i=0; i<target.size(); i++)
+      if (target.binary(i))
         if (value(i)>value_cutoff_)
           a_+=weight(i);
@@ -156 +184 @@
     // If a column sum or a row sum is zero, the table is non-sense
     if ((a_==0 || d_==0) && (c_==0 || b_==0)){
+      // Peter should throw an exception
       std::cerr << "Warning: Fisher: Table is not valid\n";
       return 1;

trunk/c++_tools/statistics/Fisher.h

@@ -118 +118 @@
     ///
     double score(const classifier::Target& target, 
+                 const classifier::DataLookupWeighted1D& value);
+
+
+    ///
+    /// Weighted version of score. Each element in 2x2 table is
+    /// calculated as \f$ \sum w_i \f$, so when each weight is
+    /// unitary the same table is created as in the unweighted version
+    ///
+    /// @return odds ratio
+    ///
+    /// @see score
+    ///
+    double score(const classifier::Target& target, 
                  const utility::vector& value,
                  const utility::vector& weight); 

trunk/c++_tools/statistics/FoldChange.cc

@@ -5 +5 @@
 #include <c++_tools/statistics/Averager.h>
 #include <c++_tools/statistics/AveragerWeighted.h>
+#include <c++_tools/classifier/DataLookupWeighted1D.h>
 #include <c++_tools/classifier/Target.h>
 
@@ -50 +51 @@
 
   double FoldChange::score(const classifier::Target& target, 
+                           const classifier::DataLookupWeighted1D& value) 
+  {
+    weighted_=true;
+    AveragerWeighted pos;
+    AveragerWeighted neg;
+
+    for (size_t i=0; i<value.size(); i++) 
+      if (target.binary(i))
+        pos.add(value.data(i),value.weight(i));
+      else
+        neg.add(value.data(i),value.weight(i));
+          
+    if (absolute_)
+      return fabs(pos.mean()-neg.mean());
+    return pos.mean()-neg.mean();
+  }
+
+
+  double FoldChange::score(const classifier::Target& target, 
                            const utility::vector& value,
                            const utility::vector& weight) 

trunk/c++_tools/statistics/FoldChange.h

@@ -30 +30 @@
     /// @return difference of the means of the two classes
     ///
-    /// @param target is +1 or -1 
+    /// @param target 
     /// @param value vector of the values
     /// 
@@ -37 +37 @@
  
     ///
+    /// @return difference of the means of the two classes
+    ///
+    /// @param target 
+    /// @param value vector of the values (with weights)
+    /// 
+    double score(const classifier::Target& target, 
+                 const classifier::DataLookupWeighted1D& value); 
+ 
+    ///
     /// @return difference of the weighted means of the two classes
     ///
-    /// @param target is +1 or -1 
+    /// @param target 
     /// @param value vector of the values
     /// @param weight vector of accompanied weight to the values

trunk/c++_tools/statistics/Pearson.cc

@@ -5 +5 @@
 #include <c++_tools/statistics/AveragerPairWeighted.h>
 #include <c++_tools/utility/vector.h>
+#include <c++_tools/classifier/DataLookupWeighted1D.h>
 #include <c++_tools/classifier/Target.h>
 
@@ -59 +60 @@
     
   double Pearson::score(const classifier::Target& target, 
+                        const classifier::DataLookupWeighted1D& value)
+  {
+    weighted_=true;
+    AveragerPairWeighted ap;
+    for (size_t i=0; i<target.size(); i++){
+      if (target.binary(i))
+        ap.add(1, value.data(i),1,value.weight(i));
+      else
+        ap.add(-1, value.data(i),1,value.weight(i));
+      nof_samples_ = target.size();
+    }
+    r_ = ap.correlation();
+    if (r_<0 && absolute_)
+      return -r_;
+     
+    return r_;
+  }
+
+  double Pearson::score(const classifier::Target& target, 
                         const utility::vector& value,
                         const utility::vector& weight)

trunk/c++_tools/statistics/Pearson.h

@@ -53 +53 @@
     ///
     double score(const classifier::Target& target, 
+                 const classifier::DataLookupWeighted1D& value); 
+
+    ///
+    /// \f$ \frac{\vert \sum_iw^2_i(x_i-\bar{x})(y_i-\bar{y})\vert }
+    /// {\sqrt{\sum_iw^2_i(x_i-\bar{x})^2\sum_iw^2_i(y_i-\bar{y})^2}}
+    /// \f$, where \f$ m_x = \frac{\sum w_ix_i}{\sum w_i} \f$ and \f$
+    /// m_x = \frac{\sum w_ix_i}{\sum w_i} \f$. This expression is
+    /// chosen to get a correlation equal to unity when \a x and \a y
+    /// are equal. @return absolute value of weighted version of
+    /// Pearson correlation.
+    ///
+    double score(const classifier::Target& target, 
                  const utility::vector& value,
                  const utility::vector& weight); 

trunk/c++_tools/statistics/ROC.cc

@@ -1 +1 @@
 // $Id$
 
-#include <c++_tools/statistics/ROC.h>
-#include <c++_tools/utility/stl_utility.h>
-#include <c++_tools/utility/vector.h>
+#include "c++_tools/statistics/ROC.h"
+
+#include "c++_tools/classifier/DataLookupWeighted1D.h"
+#include "c++_tools/utility/stl_utility.h"
+#include "c++_tools/utility/vector.h"
 
 #include <gsl/gsl_cdf.h>
@@ -97 +99 @@
     return area_;
   }
+
+
+  // Peter, should be possible to do this in NlogN
+  double ROC::score(const classifier::Target& target, 
+                    const classifier::DataLookupWeighted1D& value)
+  {
+    weighted_=true;
+
+    vec_pair_.clear();
+    vec_pair_.reserve(target.size());
+    for (unsigned int i=0; i<target.size(); i++)
+      if (value.weight(i))
+        vec_pair_.push_back(std::make_pair(target.binary(i),value.data(i)));
+
+    std::sort(vec_pair_.begin(),vec_pair_.end(),
+              utility::pair_value_compare<int, double>());
+
+    area_=0;
+    nof_pos_=0;
+    double max_area=0;
+
+    for (size_t i=0; i<n(); i++)
+      if (target.binary(i))
+        for (size_t j=0; j<n(); j++)
+          if (!target.binary(j)){
+            if (value.data(i)>value.data(j))
+              area_+=value.weight(i)*value.weight(j);
+            max_area+=value.weight(i)*value.weight(j);
+          }
+    
+    area_/=max_area;
+    
+    if (area_<0.5 && absolute_)
+      area_=1.0-area_;
+    
+    return area_;
+  }
+
 
   // Peter, should be possible to do this in NlogN

trunk/c++_tools/statistics/ROC.h

@@ -47 +47 @@
                  const utility::vector& value); 
     
+    /// Function taking values, target, weight and a vector defining
+    /// what samples to use. The area is defines as \f$ \frac{\sum
+    /// w^+w^-}{\sum w^+w^-}\f$, where the sum in the numerator goes
+    /// over all pairs where value+ is larger than value-. The
+    /// denominator goes over all pairs. If target is equal to 1,
+    /// sample belonges to class + otherwise sample belongs to class
+    /// -. @return wheighted version of area under the ROC curve. If
+    /// the area is less than 0.5 and absolute=true, 1-area is
+    /// returned. Complexity is \f$ N^2 \f$ where \f$ N \f$ is number
+    /// of samples.
+    ///
+    double score(const classifier::Target& target, 
+                 const classifier::DataLookupWeighted1D& value); 
+       
+
     /// Function taking values, target, weight and a vector defining
     /// what samples to use. The area is defines as \f$ \frac{\sum

trunk/c++_tools/statistics/SNR.cc

@@ -1 +1 @@
 // $Id$
 
-
-// Thep C++ Tools
-#include <c++_tools/statistics/SNR.h>
-#include <c++_tools/statistics/Averager.h>
-#include <c++_tools/statistics/AveragerWeighted.h>
-#include <c++_tools/classifier/Target.h>
+#include "c++_tools/statistics/SNR.h"
+#include "c++_tools/statistics/Averager.h"
+#include "c++_tools/statistics/AveragerWeighted.h"
+#include "c++_tools/classifier/DataLookupWeighted1D.h"
+#include "c++_tools/classifier/Target.h"
 
 namespace theplu {
@@ -39 +38 @@
 
   double SNR::score(const classifier::Target& target, 
+                    const classifier::DataLookupWeighted1D& value)
+  {
+    weighted_=true;
+    statistics::AveragerWeighted positive;
+    statistics::AveragerWeighted negative;
+    for(size_t i=0; i<target.size(); i++){
+      if (target.binary(i))
+        positive.add(value.data(i),value.weight(i));
+      else
+        negative.add(value.data(i),value.weight(i));
+    }
+    double diff = positive.mean() - negative.mean();
+    double denom=positive.std()+negative.std();
+    assert(denom);
+    score_=diff/denom;
+    if(positive.sum_w()==0 || negative.sum_w()==0) 
+      score_=0;
+    if (score_<0 && absolute_)
+      score_=-score_;    
+    return score_;
+  }
+
+
+
+  double SNR::score(const classifier::Target& target, 
                     const utility::vector& value,
                     const utility::vector& weight)

trunk/c++_tools/statistics/SNR.h

@@ -11 +11 @@
   namespace utility {
     class vector;
+  }
+  namespace classifier {
+    class DataLookWeighted1D;
   }
 namespace statistics {  
@@ -45 +48 @@
     ///
     double score(const classifier::Target& target, 
+                 const classifier::DataLookupWeighted1D& value);          
+
+    ///
+    /// Weighted version of SNR @return t-score if absolute=true
+    /// absolute value of t-score is returned.
+    ///
+    double score(const classifier::Target& target, 
                  const utility::vector& value, 
                  const utility::vector& weight);          

trunk/c++_tools/statistics/Score.h

@@ -17 +17 @@
   class Target;
   class DataLookup1D;
+  class DataLookupWeighted1D;
 }
 
@@ -75 +76 @@
  
     ///
+    /// Function calculating the score in a weighted fashion. In
+    /// absolute mode, also the score using negated class labels is
+    /// calculated, and the largest of the two scores are
+    /// calculated. Absolute mode should be used when two-tailed test
+    /// is wanted.
+    ///
+    virtual double 
+    score(const classifier::Target& target, 
+          const classifier::DataLookupWeighted1D& value) = 0; 
+ 
+    ///
     /// Function calculating the weighted version of score. In
     /// absolute mode, also the score using negated class labels is

trunk/c++_tools/statistics/WilcoxonFoldChange.cc

@@ -38 +38 @@
   }
 
+
+  double WilcoxonFoldChange::score(const classifier::Target& target, 
+                                   const classifier::DataLookupWeighted1D& value)
+  {
+    std::cerr << " WilcoxonFoldChange::score  not implemented" << std::endl;
+    return 0;
+  }
+
+
   double WilcoxonFoldChange::score(const classifier::Target& target, 
                                    const utility::vector& value,

trunk/c++_tools/statistics/WilcoxonFoldChange.h

@@ -38 +38 @@
     /// @return difference of the weighted means of the two classes
     ///
+    /// @param value vector of the values (with weights)
+    /// @train_set defining which values to use (number of values used
+    /// in the calculation is equal to size of \a train_set)
+    ///
+    /// @note not implemented
+    ///
+    double score(const classifier::Target& target, 
+                 const classifier::DataLookupWeighted1D& value); 
+ 
+    /// 
+    /// @return difference of the weighted means of the two classes
+    ///
     /// @param value vector of the values
     /// @param weight vector of accompanied weight to the values

trunk/c++_tools/statistics/tScore.cc

@@ -1 +1 @@
 // $Id$
 
-// System includes
+#include "c++_tools/statistics/tScore.h"
+#include "c++_tools/statistics/Averager.h"
+#include "c++_tools/statistics/AveragerWeighted.h"
+#include "c++_tools/classifier/DataLookupWeighted1D.h"
+#include "c++_tools/classifier/Target.h"
+
 #include <cassert>
 #include <cmath>
 
-// Thep C++ Tools
-#include <c++_tools/statistics/tScore.h>
-#include <c++_tools/statistics/Averager.h>
-#include <c++_tools/statistics/AveragerWeighted.h>
-#include <c++_tools/classifier/Target.h>
 
 namespace theplu {
@@ -41 +41 @@
     return t_;
   }
+
+
+  double tScore::score(const classifier::Target& target, 
+                       const classifier::DataLookupWeighted1D& value)
+  {
+    weighted_=true;
+
+    statistics::AveragerWeighted positive;
+    statistics::AveragerWeighted negative;
+    for(size_t i=0; i<target.size(); i++){
+      if (target.binary(i))
+        positive.add(value.data(i),value.weight(i));
+      else
+        negative.add(value.data(i),value.weight(i));
+    }
+    double diff = positive.mean() - negative.mean();
+    dof_=positive.n()+negative.n()-2;
+    double s2=(positive.sum_xx_centered()+negative.sum_xx_centered())/dof_;
+    t_=diff/sqrt(s2/positive.n()+s2/(negative.n()));
+    if (t_<0 && absolute_)
+      t_=-t_;
+
+    if(positive.sum_w()==0 || negative.sum_w()==0)
+      t_=0;
+    return t_;
+  }
+
 
   double tScore::score(const classifier::Target& target, 

trunk/c++_tools/statistics/tScore.h

@@ -62 +62 @@
     ///
     double score(const classifier::Target& target, 
+                 const classifier::DataLookupWeighted1D& value); 
+
+    ///
+    /// Calculates the weighted t-score, i.e. the ratio between
+    /// difference in mean and standard deviation of this
+    /// difference. \f$ t = \frac{ m_x - m_y }{
+    /// \frac{s2}{n_x}+\frac{s2}{n_y}} \f$ where \f$ m \f$ is the
+    /// weighted mean, n is the weighted version of number of data
+    /// points and \f$ s2 \f$ is an estimation of the variance \f$ s^2
+    /// = \frac{ \sum_i w_i(x_i-m_x)^2 + \sum_i w_i(y_i-m_y)^2 }{ n_x
+    /// + n_y - 2 } \f$. See AveragerWeighted for details.
+    ///
+    /// @return t-score if absolute=true absolute value of t-score
+    /// is returned
+    ///
+    double score(const classifier::Target& target, 
                  const utility::vector& value, 
                  const utility::vector& weight);