source: branches/DatasetFeatureCorrelation/HeuristicLab.Problems.DataAnalysis/3.4/OnlineCalculators/SpearmansRankCorrelationCoefficientCalculator.cs @ 8276

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2012 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 *
 * This file is part of HeuristicLab.
 *
 * HeuristicLab is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * HeuristicLab is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */
#endregion

using System;
using System.Collections.Generic;
using System.Linq;

namespace HeuristicLab.Problems.DataAnalysis {
  public class SpearmansRankCorrelationCoefficientCalculator {

    public static double Calculate(IEnumerable<double> originalValues, IEnumerable<double> estimatedValues, out OnlineCalculatorError errorState) {
      double rs = Math.Abs(Spear(originalValues, estimatedValues, out errorState));
      if (errorState != OnlineCalculatorError.None) return double.NaN;
      return rs;
    }

    /// <summary>
    /// Calculates Spearmans Rank Correlation Coefficient. Source: Numerical Recipes in C.
    /// </summary>
    private static double Spear(IEnumerable<double> xs, IEnumerable<double> ys, out OnlineCalculatorError errorState) {
      double[] xsArr = xs.ToArray();
      double[] ysArr = ys.ToArray();
      if (xsArr.Length != ysArr.Length) throw new ArgumentException("The number of elements in xs and ys does not match");

      int n = xsArr.Length;
      Array.Sort(xsArr, ysArr);
      double sf = CRank(xsArr);
      Array.Sort(ysArr, xsArr);
      double sg = CRank(ysArr);

      double d = 0.0;
      for (int j = 0; j < n; j++) //Sum the squared difference of ranks.
        d += Math.Pow(xsArr[j] - ysArr[j], 2);

      double en, en3n, aved, fac, rs; // vard, zd, t, df, probd, probrs;
      en = n;
      en3n = en * en * en - en;
      aved = en3n / 6.0 - (sf + sg) / 12.0;
      fac = (1.0 - sf / en3n) * (1.0 - sg / en3n);
      //vard = ((en - 1.0) * en * en * Math.Pow(en + 1.0, 2) / 36.0) * fac;
      //zd = (d - aved) / Math.Sqrt(vard);
      //probd = erfcc(Math.Abs(zd) / 1.4142136);
      rs = (1.0 - (6.0 / en3n) * (d + (sf + sg) / 12.0)) / Math.Sqrt(fac);
      //fac = (rs + 1.0) * (1.0 - rs);
      //if (fac > 0.0) {
      //  t = rs * Math.Sqrt((en - 2.0) / fac);
      //  df = en - 2.0;
      //  probrs = betai(0.5 * df, 0.5, df / (df + t * t));
      //} else {
      //  probrs = 0.0;
      //}
      errorState = OnlineCalculatorError.None;
      return rs;
    }

    /// <summary>
    /// Calculates midranks. Source: Numerical Recipes in C.
    /// </summary>
    /// <param name="w">Sorted array of elements, replaces the elements by their rank, including midranking of ties</param>
    /// <returns></returns>
    private static double CRank(double[] w) {
      int i = 0;
      int n = w.Length;
      double s = 0.0;
      double t;
      while (i < n - 1) {
        if (w[i + 1] > w[i]) {    // w[i+1] must be larger or equal w[i] as w must be sorted
          // not a tie
          w[i] = i + 1;
          i++;
        } else {
          int j;
          for (j = i + 1; j < n && w[j] <= w[i]; j++) ; // how far does it go (<= effectively means == as w must be sorted)
          double rank = 0.5 * (i + j - 1);
          int k;
          for (k = i; k <= j - 1; k++) w[k] = rank; // set the rank for all tied entries
          t = j - i;
          s += t * t * t - t;
          i = j;
        }
      }

      if (i == n - 1) w[n - 1] = n - 1;   // if the last element was not tied, this is its rank
      return s;
    }
  }
}