Changeset 9363 for branches/OaaS/HeuristicLab.Problems.DataAnalysis/3.4/Implementation/Classification/ThresholdCalculators

Timestamp:

04/16/13 13:13:41 (12 years ago)

Author:

spimming

Message:

#1888:

• Merged revisions from trunk

Location:

branches/OaaS

Files:

• . (modified) (2 props)
• HeuristicLab.Problems.DataAnalysis (modified) (1 prop)
• HeuristicLab.Problems.DataAnalysis/3.4/Implementation/Classification/ThresholdCalculators/AccuracyMaximizationThresholdCalculator.cs (modified) (3 diffs)
• HeuristicLab.Problems.DataAnalysis/3.4/Implementation/Classification/ThresholdCalculators/NormalDistributionCutPointsThresholdCalculator.cs (modified) (4 diffs)

branches/OaaS

@@ -21 +21 @@
 protoc.exe
 _ReSharper.HeuristicLab 3.3 Tests
+Google.ProtocolBuffers-2.4.1.473.dll
 packages

@@ -21 +21 @@
 protoc.exe
 _ReSharper.HeuristicLab 3.3 Tests
+Google.ProtocolBuffers-2.4.1.473.dll
 packages

branches/OaaS/HeuristicLab.Problems.DataAnalysis/3.4/Implementation/Classification/ThresholdCalculators/AccuracyMaximizationThresholdCalculator.cs

@@ -53 +53 @@
 
     public static void CalculateThresholds(IClassificationProblemData problemData, IEnumerable<double> estimatedValues, IEnumerable<double> targetClassValues, out double[] classValues, out double[] thresholds) {
-      int slices = 100;
-      double minThresholdInc = 10e-5; // necessary to prevent infinite loop when maxEstimated - minEstimated is effectively zero (constant model)
+      const int slices = 100;
+      const double minThresholdInc = 10e-5; // necessary to prevent infinite loop when maxEstimated - minEstimated is effectively zero (constant model)
       List<double> estimatedValuesList = estimatedValues.ToList();
       double maxEstimatedValue = estimatedValuesList.Max();
@@ -61 +61 @@
       var estimatedAndTargetValuePairs =
         estimatedValuesList.Zip(targetClassValues, (x, y) => new { EstimatedValue = x, TargetClassValue = y })
-        .OrderBy(x => x.EstimatedValue)
-        .ToList();
+        .OrderBy(x => x.EstimatedValue).ToList();
 
-      classValues = problemData.ClassValues.OrderBy(x => x).ToArray();
+      classValues = estimatedAndTargetValuePairs.GroupBy(x => x.TargetClassValue)
+        .Select(x => new { Median = x.Select(y => y.EstimatedValue).Median(), Class = x.Key })
+        .OrderBy(x => x.Median).Select(x => x.Class).ToArray();
+
       int nClasses = classValues.Length;
       thresholds = new double[nClasses];
       thresholds[0] = double.NegativeInfinity;
-      // thresholds[thresholds.Length - 1] = double.PositiveInfinity;
 
       // incrementally calculate accuracy of all possible thresholds
@@ -85 +86 @@
             //all positives
             if (pair.TargetClassValue.IsAlmost(classValues[i - 1])) {
-              if (pair.EstimatedValue > lowerThreshold && pair.EstimatedValue < actualThreshold)
+              if (pair.EstimatedValue > lowerThreshold && pair.EstimatedValue <= actualThreshold)
                 //true positive
-                classificationScore += problemData.GetClassificationPenalty(classValues[i - 1], classValues[i - 1]);
+                classificationScore += problemData.GetClassificationPenalty(pair.TargetClassValue, pair.TargetClassValue);
               else
                 //false negative
-                classificationScore += problemData.GetClassificationPenalty(classValues[i], classValues[i - 1]);
+                classificationScore += problemData.GetClassificationPenalty(pair.TargetClassValue, classValues[i]);
             }
               //all negatives
             else {
-              if (pair.EstimatedValue > lowerThreshold && pair.EstimatedValue < actualThreshold)
-                //false positive
-                classificationScore += problemData.GetClassificationPenalty(classValues[i - 1], classValues[i]);
-              else
-                //true negative, consider only upper class
-                classificationScore += problemData.GetClassificationPenalty(classValues[i], classValues[i]);
+              //false positive
+              if (pair.EstimatedValue > lowerThreshold && pair.EstimatedValue <= actualThreshold)
+                classificationScore += problemData.GetClassificationPenalty(pair.TargetClassValue, classValues[i - 1]);
+              else if (pair.EstimatedValue <= lowerThreshold)
+                classificationScore += problemData.GetClassificationPenalty(pair.TargetClassValue, classValues[i - 2]);
+              else if (pair.EstimatedValue > actualThreshold) {
+                if (pair.TargetClassValue < classValues[i - 1]) //negative in wrong class, consider upper class
+                  classificationScore += problemData.GetClassificationPenalty(pair.TargetClassValue, classValues[i]);
+                else //true negative, must be optimized by the other thresholds
+                  classificationScore += problemData.GetClassificationPenalty(pair.TargetClassValue, pair.TargetClassValue);
+              }
             }
           }

branches/OaaS/HeuristicLab.Problems.DataAnalysis/3.4/Implementation/Classification/ThresholdCalculators/NormalDistributionCutPointsThresholdCalculator.cs

@@ -53 +53 @@
 
     public static void CalculateThresholds(IClassificationProblemData problemData, IEnumerable<double> estimatedValues, IEnumerable<double> targetClassValues, out double[] classValues, out double[] thresholds) {
-      double maxEstimatedValue = estimatedValues.Max();
-      double minEstimatedValue = estimatedValues.Min();
       var estimatedTargetValues = Enumerable.Zip(estimatedValues, targetClassValues, (e, t) => new { EstimatedValue = e, TargetValue = t }).ToList();
+      double estimatedValuesRange = estimatedValues.Range();
 
       Dictionary<double, double> classMean = new Dictionary<double, double>();
@@ -72 +71 @@
         }
       }
+
       double[] originalClasses = classMean.Keys.OrderBy(x => x).ToArray();
       int nClasses = originalClasses.Length;
@@ -82 +82 @@
           // calculate all thresholds
           CalculateCutPoints(classMean[class0], classStdDev[class0], classMean[class1], classStdDev[class1], out x1, out x2);
-          if (!thresholdList.Any(x => x.IsAlmost(x1))) thresholdList.Add(x1);
-          if (!thresholdList.Any(x => x.IsAlmost(x2))) thresholdList.Add(x2);
+
+          // if the two cut points are too close (for instance because the stdDev=0)
+          // then move them by 0.1% of the range of estimated values
+          if (x1.IsAlmost(x2)) {
+            x1 -= 0.001 * estimatedValuesRange;
+            x2 += 0.001 * estimatedValuesRange;
+          }
+          if (!double.IsInfinity(x1) && !thresholdList.Any(x => x.IsAlmost(x1))) thresholdList.Add(x1);
+          if (!double.IsInfinity(x2) && !thresholdList.Any(x => x.IsAlmost(x2))) thresholdList.Add(x2);
         }
       }
       thresholdList.Sort();
+
+      // add small value and large value for the calculation of most influential class in each thresholded section
       thresholdList.Insert(0, double.NegativeInfinity);
-
-      // determine class values for each partition separated by a threshold by calculating the density of all class distributions
-      // all points in the partition are classified as the class with the maximal density in the parition
-      List<double> classValuesList = new List<double>();
-      for (int i = 0; i < thresholdList.Count; i++) {
-        double m;
-        if (double.IsNegativeInfinity(thresholdList[i])) {
-          m = thresholdList[i + 1] - 1.0; // smaller than the smalles non-infinity threshold
-        } else if (i == thresholdList.Count - 1) {
-          // last threshold
-          m = thresholdList[i] + 1.0; // larger than the last threshold
-        } else {
-          m = thresholdList[i] + (thresholdList[i + 1] - thresholdList[i]) / 2.0; // middle of partition
-        }
-
-        // determine class with maximal probability density in m
-        double maxDensity = double.MinValue;
-        double maxDensityClassValue = -1;
-        foreach (var classValue in originalClasses) {
-          double density = NormalDensity(m, classMean[classValue], classStdDev[classValue]);
+      thresholdList.Add(double.PositiveInfinity);
+
+
+      // find the most likely class for the points between thresholds m
+      List<double> filteredThresholds = new List<double>();
+      List<double> filteredClassValues = new List<double>();
+      for (int i = 0; i < thresholdList.Count - 1; i++) {
+        // determine class with maximal density mass between the thresholds
+        double maxDensity = DensityMass(thresholdList[i], thresholdList[i + 1], classMean[originalClasses[0]], classStdDev[originalClasses[0]]);
+        double maxDensityClassValue = originalClasses[0];
+        foreach (var classValue in originalClasses.Skip(1)) {
+          double density = DensityMass(thresholdList[i], thresholdList[i + 1], classMean[classValue], classStdDev[classValue]);
           if (density > maxDensity) {
             maxDensity = density;
@@ -113 +114 @@
           }
         }
-        classValuesList.Add(maxDensityClassValue);
-      }
-
-      // only keep thresholds at which the class changes 
-      // class B overrides threshold s. So only thresholds r and t are relevant and have to be kept
-      // 
-      //      A    B  C
-      //       /\  /\/\        
-      //      / r\/ /\t\       
-      //     /   /\/  \ \      
-      //    /   / /\s  \ \     
-      //  -/---/-/ -\---\-\----
-      List<double> filteredThresholds = new List<double>();
-      List<double> filteredClassValues = new List<double>();
-      filteredThresholds.Add(thresholdList[0]);
-      filteredClassValues.Add(classValuesList[0]);
-      for (int i = 0; i < classValuesList.Count - 1; i++) {
-        if (classValuesList[i] != classValuesList[i + 1]) {
-          filteredThresholds.Add(thresholdList[i + 1]);
-          filteredClassValues.Add(classValuesList[i + 1]);
-        }
-      }
+        if (maxDensity > double.NegativeInfinity &&
+          (filteredClassValues.Count == 0 || !maxDensityClassValue.IsAlmost(filteredClassValues.Last()))) {
+          filteredThresholds.Add(thresholdList[i]);
+          filteredClassValues.Add(maxDensityClassValue);
+        }
+      }
+
+      if (filteredThresholds.Count == 0 || !double.IsNegativeInfinity(filteredThresholds.First())) {
+        // this happens if there are no thresholds (distributions for all classes are exactly the same)
+        // or when the CDF up to the first threshold is zero
+        // -> all samples should be classified as the class with the most observations
+        // group observations by target class and select the class with largest count 
+        double mostFrequentClass = targetClassValues.GroupBy(c => c)
+                              .OrderBy(g => g.Count())
+                              .Last().Key;
+        filteredThresholds.Insert(0, double.NegativeInfinity);
+        filteredClassValues.Insert(0, mostFrequentClass);
+      }
+
       thresholds = filteredThresholds.ToArray();
       classValues = filteredClassValues.ToArray();
     }
 
-    private static double NormalDensity(double x, double mu, double sigma) {
+    private static double sqr2 = Math.Sqrt(2.0);
+    // returns the density function of the standard normal distribution at x
+    private static double NormalCDF(double x) {
+      return 0.5 * (1 + alglib.errorfunction(x / sqr2));
+    }
+
+    // approximation of the log of the normal cummulative distribution from the lightspeed toolbox by Tom Minka
+    // http://research.microsoft.com/en-us/um/people/minka/software/lightspeed/
+    private static double[] c = new double[] { -1, 5 / 2.0, -37 / 3.0, 353 / 4.0, -4081 / 5.0, 55205 / 6.0, -854197 / 7.0 };
+    private static double LogNormalCDF(double x) {
+      if (x >= -6.5)
+        // calculate the log directly if x is large enough
+        return Math.Log(NormalCDF(x));
+      else {
+        double z = Math.Pow(x, -2);
+        // asymptotic series for logcdf
+        double y = z * (c[0] + z * (c[1] + z * (c[2] + z * (c[3] + z * (c[4] + z * (c[5] + z * c[6]))))));
+        return y - 0.5 * Math.Log(2 * Math.PI) - 0.5 * x * x - Math.Log(-x);
+      }
+    }
+
+    // determines the value NormalCDF(mu,sigma, upper)  - NormalCDF(mu, sigma, lower)
+    // = the integral of the PDF of N(mu, sigma) in the range [lower, upper]
+    private static double DensityMass(double lower, double upper, double mu, double sigma) {
       if (sigma.IsAlmost(0.0)) {
-        if (x.IsAlmost(mu)) return 1.0; else return 0.0;
+        if (lower < mu && mu < upper) return 0.0; // all mass is between lower and upper
+        else return double.NegativeInfinity; // no mass is between lower and upper
+      }
+
+      if (lower > mu) {
+        return DensityMass(-upper, -lower, -mu, sigma);
+      }
+
+      upper = (upper - mu) / sigma;
+      lower = (lower - mu) / sigma;
+      if (double.IsNegativeInfinity(lower)) return LogNormalCDF(upper);
+
+      return LogNormalCDF(upper) + Math.Log(1 - Math.Exp(LogNormalCDF(lower) - LogNormalCDF(upper)));
+    }
+
+    // Calculates the points x1 and x2 where the distributions N(m1, s1) == N(m2,s2). 
+    // In the general case there should be two cut points. If either s1 or s2 is 0 then x1==x2. 
+    // If both s1 and s2 are zero than there are no cut points but we should return something reasonable (e.g. (m1 + m2) / 2) then.
+    private static void CalculateCutPoints(double m1, double s1, double m2, double s2, out double x1, out double x2) {
+      if (s1.IsAlmost(s2)) {
+        if (m1.IsAlmost(m2)) {
+          x1 = double.NegativeInfinity;
+          x2 = double.NegativeInfinity;
+        } else {
+          // s1==s2 and m1 != m2
+          // return something reasonable. cut point should be half way between m1 and m2
+          x1 = (m1 + m2) / 2;
+          x2 = double.NegativeInfinity;
+        }
+      } else if (s1.IsAlmost(0.0)) {
+        // when s1 is 0.0 the cut points are exactly at m1 ...
+        x1 = m1;
+        x2 = m1;
+      } else if (s2.IsAlmost(0.0)) {
+        // ... same for s2
+        x1 = m2;
+        x2 = m2;
       } else {
-        return (1.0 / Math.Sqrt(2.0 * Math.PI * sigma * sigma)) * Math.Exp(-((x - mu) * (x - mu)) / (2.0 * sigma * sigma));
-      }
-    }
-
-    private static void CalculateCutPoints(double m1, double s1, double m2, double s2, out double x1, out double x2) {
-      double a = (s1 * s1 - s2 * s2);
-      x1 = -(-m2 * s1 * s1 + m1 * s2 * s2 + Math.Sqrt(s1 * s1 * s2 * s2 * ((m1 - m2) * (m1 - m2) + 2.0 * (-s1 * s1 + s2 * s2) * Math.Log(s2 / s1)))) / a;
-      x2 = (m2 * s1 * s1 - m1 * s2 * s2 + Math.Sqrt(s1 * s1 * s2 * s2 * ((m1 - m2) * (m1 - m2) + 2.0 * (-s1 * s1 + s2 * s2) * Math.Log(s2 / s1)))) / a;
+        if (s2 < s1) {
+          // make sure s2 is the larger std.dev.
+          CalculateCutPoints(m2, s2, m1, s1, out x1, out x2);
+        } else {
+          // general case
+          // calculate the solutions x1, x2 where N(m1,s1) == N(m2,s2)
+          double g = Math.Sqrt(2 * s2 * s2 * Math.Log(s2 / s1) - 2 * s1 * s1 * Math.Log(s2 / s1) - 2 * m1 * m2 + m1 * m1 + m2 * m2);
+          double s = (s1 * s1 - s2 * s2);
+          x1 = (m2 * s1 * s1 - m1 * s2 * s2 + s1 * s2 * g) / s;
+          x2 = -(m1 * s2 * s2 - m2 * s1 * s1 + s1 * s2 * g) / s;
+        }
+      }
     }
   }

@@ -21 +21 @@
 protoc.exe
 _ReSharper.HeuristicLab 3.3 Tests
+Google.ProtocolBuffers-2.4.1.473.dll
 packages