Changeset 17787

Timestamp:

12/02/20 00:15:24 (4 years ago)

Author:

bburlacu

Message:

#3090: Implement extension to Welford's algorithm by Schubert et al in the OnlinePearsonsRCalculator.

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 OnlineCovarianceCalculator covCalculator = new OnlineCovarianceCalculator();
-    private OnlineMeanAndVarianceCalculator sxCalculator = new OnlineMeanAndVarianceCalculator();
-    private OnlineMeanAndVarianceCalculator syCalculator = new OnlineMeanAndVarianceCalculator();
+    private double sumX;
+    private double sumY;
+    private double sumWe;
+
+    private double sumXX;
+    private double sumYY;
+    private double sumXY;
+
+    private OnlineCalculatorError errorState;
 
     public double R {
       get {
-        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;
+        if (!(sumXX > 0.0 && sumYY > 0.0)) {
+          return (sumXX == sumYY) ? 1.0 : 0.0;
         }
+        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() { }
@@ -49 +100 @@
     protected OnlinePearsonsRCalculator(OnlinePearsonsRCalculator original, Cloner cloner)
       : base(original, cloner) {
-      covCalculator = cloner.Clone(original.covCalculator);
-      sxCalculator = cloner.Clone(original.sxCalculator);
-      syCalculator = cloner.Clone(original.syCalculator);
+      sumX = original.sumX;
+      sumY = original.sumY;
+      sumXX = original.sumXX;
+      sumYY = original.sumYY;
+      sumXY = original.sumXY;
+      sumWe = original.sumWe;
+      errorState = original.ErrorState;
     }
     public override IDeepCloneable Clone(Cloner cloner) {
@@ -59 +114 @@
     #region IOnlineCalculator Members
     public OnlineCalculatorError ErrorState {
-      get { return covCalculator.ErrorState | sxCalculator.PopulationVarianceErrorState | syCalculator.PopulationVarianceErrorState; }
+      get { return errorState; }
     }
     public double Value {
@@ -65 +120 @@
     }
     public void Reset() {
-      covCalculator.Reset();
-      sxCalculator.Reset();
-      syCalculator.Reset();
+      sumXX = sumYY = sumXY = sumX = sumY = sumWe = 0.0;
+      errorState = OnlineCalculatorError.None;
     }
 
     public void Add(double x, double 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);
+      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;
     }
 
@@ -80 +150 @@
 
     public static double Calculate(IEnumerable<double> first, IEnumerable<double> second, out OnlineCalculatorError errorState) {
-      IEnumerator<double> firstEnumerator = first.GetEnumerator();
-      IEnumerator<double> secondEnumerator = second.GetEnumerator();
-      var calculator = new OnlinePearsonsRCalculator();
+      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;
 
-      // 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;
+      // 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;
       }
 
-      // 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;
-      }
+      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);
     }
   }