Changeset 8658 for trunk/sources/HeuristicLab.Problems.DataAnalysis/3.4/Implementation/Classification/ThresholdCalculators/NormalDistributionCutPointsThresholdCalculator.cs

Timestamp:

09/14/12 18:48:07 (12 years ago)

Author:

gkronber

Message:

#1925 fixed some bugs in NormalDistributionCutPointsThresholdCalculator (probably more remaining)

File:

• trunk/sources/HeuristicLab.Problems.DataAnalysis/3.4/Implementation/Classification/ThresholdCalculators/NormalDistributionCutPointsThresholdCalculator.cs (modified) (5 diffs)

trunk/sources/HeuristicLab.Problems.DataAnalysis/3.4/Implementation/Classification/ThresholdCalculators/NormalDistributionCutPointsThresholdCalculator.cs

@@ -93 +93 @@
       }
       thresholdList.Sort();
-      thresholdList.Insert(0, double.NegativeInfinity);
+
+      // add small value and large value for the calculation of most influential class in each thresholded section
+      thresholdList.Insert(0, estimatedValues.Min() - 1);
+      thresholdList.Add(estimatedValues.Max() + 1);
 
       // 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>();
-      if (thresholdList.Count == 1) {
+      if (thresholdList.Count == 2) {
         // this happens if there are no thresholds (distributions for all classes are exactly the same)
         // -> all samples should be classified as the first class
@@ -105 +108 @@
         // at least one reasonable threshold ...
         // find the most likely class for the points between thresholds m
-        for (int i = 0; i < thresholdList.Count; i++) {
-          double m;
-          if (double.IsNegativeInfinity(thresholdList[i])) {
-            m = thresholdList[i + 1] - 1.0; // smaller than the smallest 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 = LogNormalDensity(m, classMean[originalClasses[0]], classStdDev[originalClasses[0]]);
+        for (int i = 0; i < thresholdList.Count - 1; i++) {
+
+          // determine class with maximal density mass between the thresholds
+          double maxDensity = LogNormalDensityMass(thresholdList[i], thresholdList[i + 1], classMean[originalClasses[0]], classStdDev[originalClasses[0]]);
           double maxDensityClassValue = originalClasses[0];
           foreach (var classValue in originalClasses.Skip(1)) {
-            double density = LogNormalDensity(m, classMean[classValue], classStdDev[classValue]);
+            double density = LogNormalDensityMass(thresholdList[i], thresholdList[i + 1], classMean[classValue], classStdDev[classValue]);
             if (density > maxDensity) {
               maxDensity = density;
@@ -139 +133 @@
       //    /   / /\s  \ \     
       //  -/---/-/ -\---\-\----
+
       List<double> filteredThresholds = new List<double>();
       List<double> filteredClassValues = new List<double>();
-      filteredThresholds.Add(thresholdList[0]);
+      filteredThresholds.Add(double.NegativeInfinity); // the smallest possible threshold for the first class
       filteredClassValues.Add(classValuesList[0]);
+      // do not include the last threshold which was just needed for the previous step
       for (int i = 0; i < classValuesList.Count - 1; i++) {
         if (!classValuesList[i].IsAlmost(classValuesList[i + 1])) {
@@ -153 +149 @@
     }
 
+    private static double LogNormalDensityMass(double lower, double upper, double mu, double sigma) {
+      if (sigma.IsAlmost(0.0)) {
+        if (lower < mu && mu < upper) return double.PositiveInfinity; // log(1)
+        else return double.NegativeInfinity; // log(0)
+      }
+
+      Func<double, double> f = (x) =>
+        x * -0.5 * Math.Log(2.0 * Math.PI * sigma * sigma) - Math.Pow(x - mu, 3) / (3 * 2.0 * sigma * sigma);
+
+      if (double.IsNegativeInfinity(lower)) return f(upper);
+      else return f(upper) - f(lower);
+    }
+
     private static double LogNormalDensity(double x, double mu, double sigma) {
-      if (sigma.IsAlmost(0.0))
-        if (mu.IsAlmost(x)) return 0.0; // (log(1))
+      if (sigma.IsAlmost(0.0)) {
+        if (x.IsAlmost(mu)) return 0.0; // log(1);
         else return double.NegativeInfinity;
-      return -0.5 * Math.Log(2.0 * Math.PI * sigma * sigma) - ((x - mu) * (x - mu)) / (2.0 * sigma * sigma);
+      }
+
+      return -0.5 * Math.Log(2.0 * Math.PI * sigma * sigma) - Math.Pow(x - mu, 2) / (2.0 * sigma * sigma);
     }
 
@@ -176 +187 @@
         x2 = m2;
       } else {
-        // scale s1 and s2 for numeric stability
-        s2 = s2 / s1;
-        s1 = 1.0;
-        double a = (s1 + s2) * (s1 - s2);
-        double g = Math.Sqrt(s1 * s1 * s2 * s2 * ((m1 - m2) * (m1 - m2) + 2.0 * (s1 * s1 + s2 * s2) * Math.Log(s2 / s1)));
-        double m1s2 = m1 * s2 * s2;
-        double m2s1 = m2 * s1 * s1;
-        x1 = -(-m2s1 + m1s2 + g) / a;
-        x2 = (m2s1 - m1s2 + g) / a;
+        if (s2 < s1) {
+          // make sure s2 is the larger std.dev.
+          CalculateCutPoints(m2, s2, m1, s1, out x1, out x2);
+        } else {
+          // scale s1 and s2 for numeric stability
+          //s2 = s2 / s1;
+          //s1 = 1.0;
+          double a = (s1 + s2) * (s1 - s2);
+          double g = Math.Sqrt(s1 * s1 * s2 * s2 * ((m1 - m2) * (m1 - m2) + 2.0 * (s1 * s1 + s2 * s2) * Math.Log(s2 / s1)));
+          double m1s2 = m1 * s2 * s2;
+          double m2s1 = m2 * s1 * s1;
+          x1 = -(-m2s1 + m1s2 + g) / a;
+          x2 = (m2s1 - m1s2 + g) / a;
+        }
       }
     }