Changeset 17788 for trunk/HeuristicLab.Problems.DataAnalysis/3.4/OnlineCalculators/OnlinePearsonsRCalculator.cs

Timestamp:

12/02/20 15:45:54 (3 years ago)

Author:

bburlacu

Message:

#3090: Revert r17787 due to numerical precision issues.

File:

• trunk/HeuristicLab.Problems.DataAnalysis/3.4/OnlineCalculators/OnlinePearsonsRCalculator.cs (modified) (5 diffs)

trunk/HeuristicLab.Problems.DataAnalysis/3.4/OnlineCalculators/OnlinePearsonsRCalculator.cs

@@ -26 +26 @@
 namespace HeuristicLab.Problems.DataAnalysis {
   public class OnlinePearsonsRCalculator : DeepCloneable, IOnlineCalculator {
-    private double sumX;
-    private double sumY;
-    private double sumWe;
-
-    private double sumXX;
-    private double sumYY;
-    private double sumXY;
-
-    private OnlineCalculatorError errorState;
+    private OnlineCovarianceCalculator covCalculator = new OnlineCovarianceCalculator();
+    private OnlineMeanAndVarianceCalculator sxCalculator = new OnlineMeanAndVarianceCalculator();
+    private OnlineMeanAndVarianceCalculator syCalculator = new OnlineMeanAndVarianceCalculator();
 
     public double R {
       get {
-        if (!(sumXX > 0.0 && sumYY > 0.0)) {
-          return (sumXX == sumYY) ? 1.0 : 0.0;
+        double xVar = sxCalculator.PopulationVariance;
+        double yVar = syCalculator.PopulationVariance;
+        if (xVar.IsAlmost(0.0) || yVar.IsAlmost(0.0)) {
+          return 0.0;
+        } else {
+          var r = covCalculator.Covariance / (Math.Sqrt(xVar) * Math.Sqrt(yVar));
+          if (r < -1.0) r = -1.0;
+          else if (r > 1.0) r = 1.0;
+          return r;
         }
-        return sumXY / Math.Sqrt(sumXX * sumYY);
       }
     }
-
-    public double MeanX { get { return sumX / sumWe; } }
-
-    public double MeanY { get { return sumY / sumWe; } }
-
-    public double NaiveCovariance { get { return sumXY / sumWe; } }
-
-    public double SampleCovariance {
-      get {
-        if (sumWe > 1.0) {
-          errorState = OnlineCalculatorError.None;
-          return sumXY / (sumWe - 1);
-        }
-        errorState = OnlineCalculatorError.InsufficientElementsAdded;
-        return double.NaN;
-      }
-    }
-
-    public double NaiveVarianceX { get { return sumXX / sumWe; } }
-
-    public double SampleVarianceX {
-      get {
-        if (sumWe > 1.0) {
-          errorState = OnlineCalculatorError.None;
-          return sumXX / (sumWe - 1);
-        }
-        errorState = OnlineCalculatorError.InsufficientElementsAdded;
-        return double.NaN;
-      }
-    }
-
-    public double NaiveStdevX { get { return Math.Sqrt(NaiveVarianceY); } }
-
-    public double SampleStdevX { get { return Math.Sqrt(SampleVarianceX); } }
-
-    public double NaiveVarianceY { get { return sumYY / sumWe; } }
-
-    public double SampleVarianceY {
-      get {
-        if (sumWe > 1.0) {
-          errorState = OnlineCalculatorError.None;
-          return sumYY / (sumWe - 1);
-        }
-        errorState = OnlineCalculatorError.InsufficientElementsAdded;
-        return double.NaN;
-      }
-    }
-
-    public double NaiveStdevY { get { return Math.Sqrt(NaiveVarianceY); } }
-
-    public double SampleStdevY { get { return Math.Sqrt(SampleVarianceX); } }
 
     public OnlinePearsonsRCalculator() { }
@@ -100 +49 @@
     protected OnlinePearsonsRCalculator(OnlinePearsonsRCalculator original, Cloner cloner)
       : base(original, cloner) {
-      sumX = original.sumX;
-      sumY = original.sumY;
-      sumXX = original.sumXX;
-      sumYY = original.sumYY;
-      sumXY = original.sumXY;
-      sumWe = original.sumWe;
-      errorState = original.ErrorState;
+      covCalculator = cloner.Clone(original.covCalculator);
+      sxCalculator = cloner.Clone(original.sxCalculator);
+      syCalculator = cloner.Clone(original.syCalculator);
     }
     public override IDeepCloneable Clone(Cloner cloner) {
@@ -114 +59 @@
     #region IOnlineCalculator Members
     public OnlineCalculatorError ErrorState {
-      get { return errorState; }
+      get { return covCalculator.ErrorState | sxCalculator.PopulationVarianceErrorState | syCalculator.PopulationVarianceErrorState; }
     }
     public double Value {
@@ -120 +65 @@
     }
     public void Reset() {
-      sumXX = sumYY = sumXY = sumX = sumY = sumWe = 0.0;
-      errorState = OnlineCalculatorError.None;
+      covCalculator.Reset();
+      sxCalculator.Reset();
+      syCalculator.Reset();
     }
 
     public void Add(double x, double y) {
-      if (sumWe <= 0.0) {
-        sumX = x;
-        sumY = y;
-        sumWe = 1;
-        return;
-      }
-      // Delta to previous mean
-      double deltaX = x * sumWe - sumX, deltaY = y * sumWe - sumY;
-      double oldWe = sumWe;
-      // Incremental update
-      sumWe += 1;
-      double f = 1.0 / (sumWe * oldWe);
-      // Update
-      sumXX += f * deltaX * deltaX;
-      sumYY += f * deltaY * deltaY;
-      // should equal weight * deltaY * neltaX!
-      sumXY += f * deltaX * deltaY;
-      // Update means
-      sumX += x;
-      sumY += y;
+      // no need to check validity of values explicitly here as it is checked in all three evaluators 
+      covCalculator.Add(x, y);
+      sxCalculator.Add(x);
+      syCalculator.Add(y);
     }
 
@@ -150 +80 @@
 
     public static double Calculate(IEnumerable<double> first, IEnumerable<double> second, out OnlineCalculatorError errorState) {
-      var x = first.GetEnumerator(); x.MoveNext();
-      var y = second.GetEnumerator(); y.MoveNext();
-      double sumXX = 0.0, sumYY = 0.0, sumXY = 0.0;
-      double sumX = x.Current, sumY = y.Current;
-      int i = 1;
+      IEnumerator<double> firstEnumerator = first.GetEnumerator();
+      IEnumerator<double> secondEnumerator = second.GetEnumerator();
+      var calculator = new OnlinePearsonsRCalculator();
 
-      // Inlined computation of Pearson correlation, to avoid allocating objects
-      // This is a numerically stabilized version, avoiding sum-of-squares.
-      while (x.MoveNext() & y.MoveNext()) {
-        double xv = x.Current, yv = y.Current;
-        // Delta to previous mean
-        double deltaX = xv * i - sumX, deltaY = yv * i - sumY;
-        // Increment count first
-        double oldi = i; // Convert to double!
-        ++i;
-        double f = 1.0 / (i * oldi);
-        // Update
-        sumXX += f * deltaX * deltaX;
-        sumYY += f * deltaY * deltaY;
-        // should equal deltaY * deltaX!
-        sumXY += f * deltaX * deltaY;
-        // Update sums
-        sumX += xv;
-        sumY += yv;
+      // always move forward both enumerators (do not use short-circuit evaluation!)
+      while (firstEnumerator.MoveNext() & secondEnumerator.MoveNext()) {
+        double original = firstEnumerator.Current;
+        double estimated = secondEnumerator.Current;
+        calculator.Add(original, estimated);
+        if (calculator.ErrorState != OnlineCalculatorError.None) break;
       }
 
-      errorState = OnlineCalculatorError.None;
-      // One or both series were constant:
-      return !(sumXX > 0.0 && sumYY > 0.0) ? sumXX == sumYY ? 1.0 : 0.0 : //
-          sumXY / Math.Sqrt(sumXX * sumYY);
+      // check if both enumerators are at the end to make sure both enumerations have the same length
+      if (calculator.ErrorState == OnlineCalculatorError.None &&
+           (secondEnumerator.MoveNext() || firstEnumerator.MoveNext())) {
+        throw new ArgumentException("Number of elements in first and second enumeration doesn't match.");
+      } else {
+        errorState = calculator.ErrorState;
+        return calculator.R;
+      }
     }
   }