Index: branches/DataAnalysis Refactoring/HeuristicLab.Algorithms.DataAnalysis/3.4/Linear/LinearDiscriminantAnalysis.cs
===================================================================
--- branches/DataAnalysis Refactoring/HeuristicLab.Algorithms.DataAnalysis/3.4/Linear/LinearDiscriminantAnalysis.cs	(revision 5664)
+++ branches/DataAnalysis Refactoring/HeuristicLab.Algorithms.DataAnalysis/3.4/Linear/LinearDiscriminantAnalysis.cs	(revision 5678)
@@ -107,9 +107,102 @@
       addition.AddSubTree(cNode);
 
-      var model = new SymbolicDiscriminantFunctionClassificationModel(tree, new SymbolicDataAnalysisExpressionTreeInterpreter(), classValues);
+
+      var model = LinearDiscriminantAnalysis.CreateDiscriminantFunctionModel(tree, new SymbolicDataAnalysisExpressionTreeInterpreter(), problemData, rows);
       SymbolicDiscriminantFunctionClassificationSolution solution = new SymbolicDiscriminantFunctionClassificationSolution(model, problemData);
+
       return solution;
     }
     #endregion
+
+    private static SymbolicDiscriminantFunctionClassificationModel CreateDiscriminantFunctionModel(ISymbolicExpressionTree tree,
+      ISymbolicDataAnalysisExpressionTreeInterpreter interpreter,
+      IClassificationProblemData problemData,
+      IEnumerable<int> rows) {
+      string targetVariable = problemData.TargetVariable;
+      List<double> originalClasses = problemData.ClassValues.ToList();
+      int nClasses = problemData.Classes;
+      List<double> estimatedValues = interpreter.GetSymbolicExpressionTreeValues(tree, problemData.Dataset, rows).ToList();
+      double maxEstimatedValue = estimatedValues.Max();
+      double minEstimatedValue = estimatedValues.Min();
+      var estimatedTargetValues =
+         (from row in problemData.TrainingIndizes
+          select new { EstimatedValue = estimatedValues[row], TargetValue = problemData.Dataset[targetVariable, row] })
+         .ToList();
+
+      Dictionary<double, double> classMean = new Dictionary<double, double>();
+      Dictionary<double, double> classStdDev = new Dictionary<double, double>();
+      // calculate moments per class
+      foreach (var classValue in originalClasses) {
+        var estimatedValuesForClass = from x in estimatedTargetValues
+                                      where x.TargetValue == classValue
+                                      select x.EstimatedValue;
+        double mean, variance;
+        OnlineMeanAndVarianceCalculator.Calculate(estimatedValuesForClass, out mean, out variance);
+        classMean[classValue] = mean;
+        classStdDev[classValue] = Math.Sqrt(variance);
+      }
+      List<double> thresholds = new List<double>();
+      for (int i = 0; i < nClasses - 1; i++) {
+        for (int j = i + 1; j < nClasses; j++) {
+          double x1, x2;
+          double class0 = originalClasses[i];
+          double class1 = originalClasses[j];
+          // calculate all thresholds
+          CalculateCutPoints(classMean[class0], classStdDev[class0], classMean[class1], classStdDev[class1], out x1, out x2);
+          if (!thresholds.Any(x => x.IsAlmost(x1))) thresholds.Add(x1);
+          if (!thresholds.Any(x => x.IsAlmost(x2))) thresholds.Add(x2);
+        }
+      }
+      thresholds.Sort();
+      thresholds.Insert(0, double.NegativeInfinity);
+      thresholds.Add(double.PositiveInfinity);
+      List<double> classValues = new List<double>();
+      for (int i = 0; i < thresholds.Count - 1; i++) {
+        double m;
+        if (double.IsNegativeInfinity(thresholds[i])) {
+          m = thresholds[i + 1] - 1.0;
+        } else if (double.IsPositiveInfinity(thresholds[i + 1])) {
+          m = thresholds[i] + 1.0;
+        } else {
+          m = thresholds[i] + (thresholds[i + 1] - thresholds[i]) / 2.0;
+        }
+
+        double maxDensity = 0;
+        double maxDensityClassValue = -1;
+        foreach (var classValue in originalClasses) {
+          double density = NormalDensity(m, classMean[classValue], classStdDev[classValue]);
+          if (density > maxDensity) {
+            maxDensity = density;
+            maxDensityClassValue = classValue;
+          }
+        }
+        classValues.Add(maxDensityClassValue);
+      }
+      List<double> filteredThresholds = new List<double>();
+      List<double> filteredClassValues = new List<double>();
+      filteredThresholds.Add(thresholds[0]);
+      filteredClassValues.Add(classValues[0]);
+      for (int i = 0; i < classValues.Count - 1; i++) {
+        if (classValues[i] != classValues[i + 1]) {
+          filteredThresholds.Add(thresholds[i + 1]);
+          filteredClassValues.Add(classValues[i + 1]);
+        }
+      }
+      filteredThresholds.Add(double.PositiveInfinity);
+
+      return new SymbolicDiscriminantFunctionClassificationModel(tree, interpreter, filteredClassValues, filteredThresholds);
+    }
+
+    private static double NormalDensity(double x, double mu, double sigma) {
+      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);
+      double b = (m1 * s2 * s2 - m2 * s1 * s1);
+      double c = 2 * s1 * s1 * s2 * s2 * Math.Log(s2) - 2 * s1 * s1 * s2 * s2 * Math.Log(s1) - s1 * s1 * m2 * m2 + s2 * s2 * m1 * m1;
+      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;
+    }
   }
 }
